Coarse-grained molecular dynamics: 
Nonlinear finite elements and finite temperature 
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Coarse-grained molecular dynamics (CGMD) is a technique developed as a concurrent multiscale 
model that couples conventional molecular dynamics (MD) to a more coarse-grained description of 
the periphery. The coarse-grained regions are modeled on a mesh in a formulation that generalizes 
conventional finite element modeling (FEM) of continuum elasticity. CGMD is derived solely from 
the MD model, however, and has no continuum parameters. As a result, it provides a coupling 
that is smooth and provides control of errors that arise at the coupling between the atomistic and 
coarse-grained regions. In this article, we elaborate on the formulation of CGMD, describing in 
detail how CGMD is applied to anharmonic solids and finite temperature simulations. As tests of 
CGMD, we present in detail the calculation of the phonon spectra for solid argon and tantalum in 
3D, demonstrating how CGMD provides a better description of the elastic waves than that provided 
by FEM. We also present elastic wave scattering calculations that show the elastic wave scattering 
is more benign in CGMD than FEM. We also discuss the dependence of scattering on the properties 
of the mesh. We introduce a rigid approximation to CGMD that eliminates internal relaxation, 
similar to the Quasicontinuum technique, and compare it to the full CGMD. 

PACS numbers: 61.43.Bn, 02.70.Ns, 62.20.-x, 62.30.-(-d 
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I. COUPLING OF LENGTH SCALES 

The Science of Materials is foremost a study of struc- 
ture. Once structure is determined other important is- 
sues such as dynamics and kinetics may be addressed. 
Structure in materials is most effectively analyzed ac- 
cording to its length scale. Materials structure at differ- 
ent scales such as crystal structure, crystal defect struc- 
ture, microstructure and macrostructure has led to the 
development of models at the atomic scale, nanoscale, mi- 
croscale, and macroscale, not to mention the mesoscale 
and a vast array of other distinctions in scale. These 
models work because the physics at one scale decouples 
to a large extent from that at other scales, provided there 
exists a sufRcient separation of scales. Then physical 
properties calculated at one scale may be passed to the 
next higher scale in a hierarchical approach that can be 
very effective 

There are systems of interest that are inherently multi- 
scale where the physics at one scale is strongly coupled to 
that at other scalesn^ Turbulence is an excellent example, 
where energy input through stirring at the macroscopic 
scale cascades down through vorticity across a range of 
length scales until it is ultimately dissipated at the short- 
est length scales. The size of the vortices varies contin- 
uously, and while there are length scales with distinct 
physics, the boundaries between them are blurred. As a 
result hierarchical models have been largely unsuccessful, 
and turbulence remains a hard problem.* This situation 
is in marked contrast to low Reynolds number flow, in 
which the physics at small length scales can be encoded 
in a few parameters, which may be computed and then 



fed into simulations at a larger scale. In this way, it has 
been possible to start with ah initio electronic structure 
calculations of H2O and through a sequence of scales end 
up with a description of tides in Buzzard's Bay:^ 

Many other examples of strongly coupled multiscale 
systems exist. Ironically, the advent of Nanoscience and 
the current focus on structures of one particular scale, the 
nanoscale, has led to the need to understand a class of 
strongly coupled multiscale systems. Consider epitaxial 
quantum dots, for example^iS The quantum dot con- 
sists of a dome of semiconductor that forms during het- 
eroepitaxy. The dome itself is typically a few nanome- 
ters to tens of nanometers across, but its size, shape 
and location are affected by the presence of other struc- 
tures during growth, even those microns away. Another 
example is a Nano-Electro-Mechanical System (NEMS) 
resonator .^iiSiiiiiS It consists of a semiconductor bar 
about 50 nm wide and a fraction of a micron long at- 
tached to the substrate at both ends. The bar itself is 
clearly nanoscale, and yet as it resonates, the oscillat- 
ing strain fields extend far into the substrate. Of course, 
there are many other examples in metallurgy and solid 
state physics, as indicated above. Remarkably, similar 
effects are beginning to be appreciated in the study of 
soft materials, chemistry and even biology. 

These systems are examples of what could be termed 
embedded nanomechanics.-^^, The mechanical properties 
of the nanoscale structure is clearly important, and these 
properties may be quite different from what would be 
predicted according to conventional macroscopic mechan- 
ics, but also important is the way the nanoscale struc- 
ture is coupled to its surroundings. Embedded nanome- 
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Resonator Simulation: 
Multiscale Domain 
Decomposition 




Coai-se Grained Molecular Dynamics \ 




FIG. 1: (color online) Schematic diagram of a coarse-grained 
simulation of a NEMS silicon microresonator.*"^ The coarse- 
grained (CG) region comprises most of the volume, but the 
molecular dynamics (MD) region contains most of the sim- 
ulated degrees of freedom. Each sphere shown in the MD 
region represents an atom. Note that the CG mesh is refined 
to the atomic scale where it joins with the MD lattice. 



chanical structures are too small to be modeled reliably 
with conventional continuum elastic theory and finite el- 
ements, and too large to be modeled by conventional 
atomistics. Even in single crystals, sub-micron dynami- 
cal regions bounded by surfaces or interfaces are affected 
by Angstrom- and nano-scale physics which causes de- 
viations from continuum elastic theory; dynamical re- 
gions larger than 0.1 micron cubed exceed the current 
limit of about one billion atoms for atomistic simulation 
of solids on a supercomputer The atomistic effects are 
compounded in materials with local defects or cracks that 
couple to long-range strain fieldsji^ The situation is not 
entirely intractable, however, because the most impor- 
tant atomistic effects are often localized to small regions 
of the system: surfaces, defects, regions of large deflection 
or internal strain, and regions of localized heating per- 
haps due to friction. The challenge is to develop a robust 
model for such an inhomogeneous system which captures 
the important atomistic effects without the prohibitive 
computational cost of a brute force atomistic simulation 
for the entire system. In this article we focus on the link 
between the micron scale and the nanoscale and develop 
a model, coarse-grained molecular dynamics (CGMD),^'' 
which bridges the disparate scales seamlessly. 

The choice to use atomistic models at the finest res- 
olution is motivated in some cases by the fact that the 
inherent length scale of the process of interest is the in- 
teratomic spacing and in other cases by the ability to 
derive interatomic potentials from quantum mechanics 



and hence built a model from first principles. Yet an- 
other motivation is that the processes of interest may be 
thermally activated, and molecular dynamics provides a 
means to simulate the thermal effects directly. Entropic 
and thermal effects are often paramount in soft matter 
systems, and in hard matter thermal activation is im- 
portant in defect diffusion, the motion of dislocations in 
metals with high Peierls barriers and many nucleation 
phenomena. Temperature is important in other ways, 
too, such as in inducing phase transitions. Also the pop- 
ulation of phonons increases with temperature, causing 
thermal expansion, changes (typically softening) in the 
clastic constants and dislocation drag at high strain rates. 
These are but a few well-known examples of the impor- 
tant role temperature plays, and thus in our development 
of multiscale models we search for methodologies capable 
of handling non-zero temperatures. 

The variation of the strain field/atomic displacements 
in inhomogeneous solid systems suggests the use of dif- 
ferent computational methodologies for different regions, 
as we mentioned above. The challenge is to meld 
them into a seamless, monolithic simulation. The first 
such proposal implemented a coupling between molecu- 
lar dynamicsiSiiS (MD) and a finite element modelSfliSi 
(FEM) implementation of continuum elastic theory us- 
ing stress consistency as the boundary condition at the 
interfaced More recently, a dynamical instability in the 
original formulation has been eliminated through the use 
of a mean force boundary condition together with uni- 
form symplectic time evolution. •^^ In both of these formu- 
lations, the same constitutive relation is used regardless 
of the size of the cells in the FEM mesh, leading to a 
discontinuity at the atomic limit. 

At its heart the FEM description of such a system re- 
lies on the ability to improve the accuracy of the sim- 
ulation by going to a finer mesh.^"'^ A mesh of varying 
coarseness is chosen, adaptively or by fiat, such that no 
single region contributes disproportionately to the error. 
These errors typically result from large strain, velocity, 
or other gradients which violate the discrete expression 
for the integral of the elastic energy density of a contin- 
uous medium. This approximation may be improved by 
mesh refinement. There is a limit, however. As the mesh 
size approaches the atomic scale, the constitutive equa- 
tions have significant errors because the expression for 
the elastic energy does not represent localized bonds and 
the standard distributed mass expression for the kinetic 
energy does not account for the fact that essentially all of 
the mass is localized in the nuclei, at least four order of 
magnitude smaller than the interatomic spacing. At this 
point, the physics of the governing equations is wrong, 
and further mesh refinement does not help. 

The approach of Refs. and|23 to improve this situa- 
tion replaces the FEM equations of motion on regions of 
the mesh that are atomic-sized with MD equations of mo- 
tion and implements a hand shaking between the MD and 
FEM regions. Although this technique is remarkably suc- 
cessful, the union is not perfectly seamless. In the FEM 
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cells approaching the atomic limit, the energy density 
varies smoothly within each cell, whereas on the other 
side of the interface, the MD energy is effectively local- 
ized to interatomic bonds. The short-wavelength modes 
of the system are able to probe this discrepancy, leading 
to errors that grow with the wavenumber. 

The quasicontinuiim technique^i^^S*^ offers another 
approach to this problem. It is a zero temperature re- 
laxation technique in which the clastic energy used in 
the FEM region is computed by applying the FEM in- 
terpolated displacement field (through the Cauchy-Born 
rule) to a reference system of atoms interacting by MD 
forces. This is a very nice idea, but it has a number of dif- 
ficulties in practice. The atoms in the reference system 
are taken to be at their zero temperature locations-no 
fluctuations are allowed. Thus, many degrees of free- 
dom are summarily set to zero (although finite temper- 
ature versions of the quasicontinuum technique are cur- 
rently under development^^). Also, the implementation 
of the quasicontinuum technique suffers from disconti- 
nuities ("ghost forces") due to the mismatch of the dis- 
placed reference systems from cell to cell and non-locality 
of atomic bondingiS^ 

Other concurrent multiscale models have been pro- 
posed recently. There are several nice review articles on 
this subject to which we direct interested readers 
The relationships of several multiscale models to CGMD 
are of particular note. A finite temperature dynamical 
model based on renormalization group ideas has been 
proposed by Curtarolo and Ceder^ The finite tempera- 
ture coarse-graining approach based on Monte Carlo cal- 
culations has been developed recently by Wu et al Ji The 
bridging scale decomposition is another approach for cou- 
pling atomistic and continuum simulations due to Wag- 
ner and Liu that provides a coupling between FEM and 
atomistics that does not require the FEM mesh to be 
refined to the atomic level^^, and it has been applied 
to finite temperature simulations by Park et al,-^ The 
projection techniques of the bridging scale decomposi- 
tion and the ensemble averages that they use are closely 
related to the techniques of CGMD introduced earlier. 
Also, we note that the assumption that no defect (dislo- 
cation) propagation takes place from an atomistic region 
into a finite element or coarse-grained region has been 
relaxed through the development of the CADD method, 
albeit so far just in two dimensions 

We have proposed a replacement for conventional fi- 
nite elements suitable for a mesh which is atomic sized in 
some regions.^'' This technique, CGMD, effectively pro- 
vides the scale-dependent constitutive equations needed 
at the interface. In the atomic limit it is guaranteed to 
reproduce the atomistic equations of motion. This en- 
ables MD regions to be coupled seamlessly to regions of 
generalized FEM, bringing the full power of MD to bear 
on important parts of the system without the compu- 
tational overhead of MD in other large, but physically 
less complex regions. The CGMD procedure is based on 
a statistical coarse graining prescription. While various 



aspects of CGMD have been introduced previously, this 
is the first article to present the model in great detail. 

This kind of multi-scale simulation poses a number of 
challenges. First, the model must have a well-behaved, 
physical response to stationary strain fields, slowly vary- 
ing in position, that extend into the CG region-there 
should be no ghost forces. Second, the model must have 
sensible thermodynamics in equilibrium. The effect of 
short- wavelength modes cannot be set to zero unless their 
energy is well above the thermal energy. Third, the sys- 
tem must have realistic dynamics, free of pathologies 
due to infiuences in the central MD region propagating 
out to unphysical interfaces, reflecting and propagating 
back into the central region to cause unphysical effects. 
Fourth, the model should exhibit well-behaved nonequi- 
librium thermodynamics, with a sensible response when 
low-lying modes are driven out of equilibrium. Finally, 
the methodology needs to be amenable to a practical im- 
plementation in terms of being able to utilize the broad 
spectrum of MD models in use, including many-body in- 
teratomic potentials that extend beyond nearest neigh- 
bors, computationally efficient domain decompositions 
for parallel distributed memory computers, visualization, 
etc. In this article we describe in detail how CGMD is 
implemented in order to meet these challenges. 

In particular, we provide the first detailed descrip- 
tion of the way CGMD may be applied to anharmonic 
solids and finite temperature simulations. Previously, we 
have described how CGMD is formulated for harmonic 
latticesii, including finite temperature contributions. We 
have also indicated how anharmonic effects are included, 
but the details have not been given. We give the details 
here, and explore their implications. We also introduce 
a rigid approximation to CGMD that eliminates internal 
relaxation, simplifying the formulation and reducing the 
computation cost of CGMD. We examine the physical 
difference between CGMD and its rigid approximation. 
As tests of CGMD, the calculation of the phonon spectra 
for solid argon and tantalum in 3D are presented. These 
calculations provide a test of the quality of the represen- 
tation of elastic wave energetics. We also present elas- 
tic wave scattering calculations, a test of how CGMD 
behaves on an irregular mesh. Whenever possible we 
present analytic formulas that contain a wealth of in- 
formation about the performance of CGMD in as much 
generality as possible (explicitly showing dependence on 
the interatomic potential, atomic masses, crystal lattice 
and mesh structure). Of course, realistic simulations in- 
volve numerical assembly and integration of the CGMD 
equations of motion. 



II. COARSE GRAINING PRESCRIPTION 

Consider a system of MD atoms in a solid, crystalline 
or amorphous, and a coarse-grained (CG) mesh partition- 
ing the solid into cells (cf. Fig. The mesh size may 
vary, so that in important regions a mesh node is assigned 
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to each equilibrium atomic position, whereas in other re- 
gions the cehs contain many atoms and the nodes need 
not coincide with atomic sites. CGMD offers a way to 
reduce the atomistic coordinates to a much smaher set of 
degrees of freedom associated with the displacement field 
at the nodes of the CG mesh, and the equations of mo- 
tion for this mean displacement field. In particular, the 
energy functional for the CG system is defined as a con- 
strained ensemble average of the atomistic energy under 
fixed thermodynamic conditions. The equations of mo- 
tion are Hamilton's equations for this conserved energy 
functional and in principle additional random and dissi- 
pative forces due to fluctuations. 

The classical ensemble must obey the constraint that 
the position and momenta of the atoms are consistent 
with the mean displacement and momentum fields. Let 
the displacement of atom /i be = — x^o where x^o 
is its equilibrium position. The displacement of mesh 
node j is a weighted average of the atomic displacements 



(1) 



where /jYi is a weighting function, related to the micro- 
scopic analog of FEM interpolating functions below. An 
analogous relation applies to the momenta . Since the 
nodal displacements are fewer or equal to the atomic po- 
sitions in number, fixing the nodal displacements and 
momenta does not (necessarily) determine the atomic 
coordinates entirely. Some subspace of phase space re- 
mains, corresponding to degrees of freedom that are miss- 
ing from the mesh. We define the GG energy as the aver- 
age energy of the canonical ensemble on this constrained 
phase space: 

£:(ufe,Ufe) = ( Hmd >ufe,u, (2) 
= J dx^dp^ Hmd e-""'"^ A/Z, (3) 



Z(uk,Uk) 



" A, 



(4) 



(5) 



where (3 = l/{kT) is the inverse temperature, Z is the 
partition function and (5(u) is a three-dimensional delta 
function. The delta functions enforce the mean field con- 
straint (Q. Note that Latin indices, j,k, . . ., denote mesh 
nodes and Greek indices, ^^v, . . ., denote atoms. The en- 
ergy Q is computed below (Eq. if^ '). 

When the mesh nodes and the atomic sites are identi- 
cal, fj^ = and the GGMD equations of motion agree 
with the atomistic equations of motion. '^^ As the mesh 
size increases some short-wavelength degrees of freedom 
are not supported by the coarse mesh. These degrees of 
freedom are not neglected entirely, because their thermo- 
dynamic average effect has been retained. This approx- 
imation is expected to be good provided the system is 



initially in thermal equilibrium, and changes to the sys- 
tem would only produce adiabatic changes in the missing 
degrees of freedom. In particular, the relaxation time of 
those degrees of freedom should be fast compared to the 
driving forces in the CG region. As long as this condition 
is satisfied, the long wavelength modes may be driven out 
of equilibrium without problems;^ 

We have written the CG energy as an internal energy, 
a function of the entropy, S, rather than the tempera- 
ture. This is designed for systems in which the short 
wavelength modes change adiabatically. This is a good 
approximation, for example, when long wavelength elas- 
tic waves propagate through a solid in the linear regime at 
finite temperaturei^^iS In other systems, the short wave- 
length modes may be in contact with a heat bath, so that 
their evolution is isothermal rather than isentropic. For 
example, the electron gas in metals can act as a heat bath 
on time scales longer than the thermal relaxation time. 
Then the Helmholz free energy. 



F(ufc,Ufe) = -A:T logZ, 



(6) 



should be used rather than the internal energy. In this 
case, the ensemble average behavior of the CG collective 
modes is exactly the same as that of the corresponding 
averaged atomic modes in the underlying atomistic sys- 
tem: 

(u.---u,J ^ /du,du,K...u,Je-/'- (7) 

j dUf^dpf, u^, • • • e-^^'^° , (8) 

which follows from plugging in the expression for F © 
and Q into {Tj) and integrating the delta functionsi^ 
This equation shows the equivalence of all unnormalized 
correlation functions, but since the partition functions 
(zero point functions) are identical, the normalized cor- 
relation functions are the same, as well. The emergence of 
the canonical distribution in other cases requires a treat- 
ment of thermal relaxation processes (cf. Section 0). It 
should be noted that even in the isothermal ensemble the 
faithfulness of correlation functions applies only to equal- 
time correlation functions at equilibrium, and consider- 
ation of dissipative processes are needed to reproduce 
interesting correlation functions such as the autocorrela- 
tion function {ui(0)ui{t)) associated with the fluctation- 
dissipation theorem. 

To end this section, we note that the definition of the 
CGMD energy may appear to neglect the well-known 
quantum mechanical contributions to lattice dynamics. 
Phonons are bosons, after all, and they should obey 
Bose-Einstein statistics. The definition of the CGMD 
energy © is clearly a classical expression based on Boltz- 
mann statistics. To what extent can it be expected to be 
valid? The reason the classical expression is valid for 
most of the conceivable applications of CGMD is that 
the Bose-Einstein distribution most strongly affects the 
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lowest states; i.e., exactly the states that are retained ex- 
plicitly in the CGMD Hamiltonian. The higher energy 
states have low occupation in equilibrium, and are not 
affected significantly by strong quantum effects such as 
Bose condensation. The CGMD Hamiltonian is there- 
fore expected to be a good description of the coarse- 
grained system. It may be necessary to use a path in- 
tegral, or other quantum mechanical version, of MD to 
treat the retained degrees of freedom at sufficiently low 
temperature^ but the internal degrees of freedom are 
described well by Eq. 



III. SHAPE FUNCTIONS 

In addition to the general framework we have presented 
for CGMD, a specific choice of the weighting functions 
is required for calculations. They result from the intro- 
duction of a set of shape functions {iVj(x)}^'°'^'' on the 
mesh from which the interpolated fields are constructed. 
The shape functions have the following properties: 



i. Wj(xfc) = 6jk, 

ii. = 1' 

J = l 

iii. C° continuity. 

The first property states that the functions arc normal- 
ized and local on the mesh nodes, x^,. The second states 
that the functions form a partition of unity, so the cen- 
ter of mass mode is represented. The third states that 
the functions are continuous but their derivatives need 
not be. This continuity guarantees that the elastic en- 
ergy, proportional to an integral of the square of the 
strain {du)sym, is well-behaved in the continuum limit. 
The interpolated displacement field is then defined by 
= A'j (x). Often there are additional consid- 

erations, such as the need to refine the mesh onto a par- 
ticular crystal lattice at the MD/CG interface^i 

Given any set of atomic displacements we can find the 
displacement field represented on the CG mesh which 
best fits these data in the least squares sense: 



E 



(9) 



where iVj^ = iVj(x^o)- This error is minimized by 



N, 



(11) 



[cf. This equation defines the weighting function /j^ 
of in terms of the interpolating function Nj{'x.). We 



note that recently this relationship introduced in CGMD 
has been generalized for use in the bridging scale, and 
other projection techniquesiSS 

The formulation we have described is appropriate to 
retain the low-lying acoustic phonon modes in the coarse- 
grained system. In some cases it is desirable to retain the 
long-wavelength optical phonons, as well. For example, 
in the study of III-V epitaxial quantum dots, internal re- 
laxation of the zinc blende structure in the strained dots 
leads to important changes in the optical spectra of the 
dotsj^^ If it is important to model the optical phonons or 
to capture the internal relaxation in a crystal lattice with 
a multiple-atom basis, each interpolation function should 
carry a band index, a, in addition to the nodal index, j: 
ivj"^ (x) . Then the basis requirements are somewhat dif- 
ferent. The functions should be local and normalized 
within each band. They should be C'' continuous apart 
from variations with the unit cell. And they should form 
a generalization of the partition of unity. In particular, 
the requirement of forming a partition of unity is the 
requirement that uniform displacement of the system be 
represented in the basis of shape functions. That transla- 
tion invariance is responsible for the k = acoustic-mode 
phonons having zero energy. We generalize the partition 
of unity requirement to the requirement that all of the 
k = phonon, both acoustic and optical, be represented. 
In particular, denote the displacement associated with 
the k = phonon for band a as u^"'', normalized such 
that 



E 



(a) 



/iGunit cell 



basis 



(12) 



where iVbasis is the number of atoms in the Wigner-Seitz 
unit cell. Then the shape functions can be defined as 



iVf'(x,,o) = u('^)iV,(x^o) 



(13) 



where Nj{x.) is a conventional shape function, such as 
linear interpolation. The generalized partition of miity 
requirement is that 



E^- 

i 



(a) 



(Xpo) 



(14) 



Note that in the case of a monatomic unit cell this shape- 
function basis is a linear combination of the shape func- 
tions we have discussed above. In that case u|f is the 
same for all lattice sites, and the orthonormal vectors 
corresponding to the three acoustic-mode phonons span 
three-dimensional space. We do not discuss an example 
of a polyatomic CGMD including optical phonons ex- 
plicitly, but the CGMD formalism continues to work. It 
should be emphasized, however, that even in polyatomic 
materials this band-index extension may not be needed 
to capture the mechanical response of interest. 



6 



IV. THE CGMD HAMILTONIAN 



We now turn to the calculation of the CGMD energy. 



A. Harmonic Lattices 

The CG energy Q may be computed in closed form 
using analytic techniques in the case of a harmonic lat- 
tice. The expression was originally presented in Ref. 0. 
We take the form of the atomistic Hamiltonian to be 



H 



MD 



D„ 



(15) 



where E'jf^ is the cohesive energy of atom /i and is 
the dynamical matrix. It acts as a tensor on the com- 
ponents of the displacement vector at each site. We re- 
express the CG energy using a parametric derivative 
of the log of the constrained partition function Q), 

E{nk,iik) = -5^1ogZ(ufc,Ufc;/3), (16) 

and we introduce the Fourier transform representation 
of the delta function (a form of Lagrange multiplier) to 
simplify the constraint (|3J) 



3Af„ 



„iA3-(uj-/j>u„) 



27r J ' 



where here, and in what follows unless stated otherwise, 
the repeated indices are summed. Expressed in this way, 
the constrained partition function for the harmonic lat- 
tice is a Gaussian integral. The complicated domain of 
integration in JSJ resulting from the constraints is re- 
placed by a simple domain plus some extra integrals. 
This standard technique gives an expression which may 
be evaluated in closed form. 

For pedagogical purposes we present the calculation 
of the CG potential energy here in some detail. The 
calculations of each of the other CG energy terms, the CG 
kinetic and anharmonic potential terms, follow the same 
basic approach. It may be helpful therefore to present 
the calculation of the CG harmonic potential energy in 
detail, and readers who are not interested in this algebra 
may skip ahead to the next paragraph at Eq. (|29|l . 

In order to get a closed-form expression for the CG 
potential energy, we make use of the well-known exact 
formula for the integral over all space of a Gaussian times 
an arbitrary polynomial prefactor. In the interest of com- 
pact notation, we combine the atomic and spatial indices 
and consider the displacements and the dynamical ma- 
trix to be objects in SA^atom-dimensional space. Similarly, 
we take all CG variables to live in 37Vnodo-dimensional 
space. The form for the Gaussian integral is known for 



these high dimensional spaces, and we evaluate the in- 
tegral |@J) through two successive Gaussian integrations: 
first an integral over the MD phase space and then an 
integral over the Lagrange multiplier space. Using l|16|l . 
we need to calculate the CG partition function. It factor- 
izes into kinetic and potential parts, Z = Zkin^pot- We 
focus on the potential energy part of the CG partition 
function: 



Zpot{uk\l3) 



(18) 



where du = ((iu)^^''*°» and d\ = (|^) Here, and 

throughout this derivation, we suppress the cohesive en- 
ergy by choosing the zero of energy such that the cohesive 
energy is zero; we then restore the cohesive energy in the 
final formulas. We first compute the integral over du by 
completing the square in the argument of the exponen- 
tial. Let 



(19) 



where we assume that the matrix inverse exists after 
a suitable regularization to deal with the zero eigenvalues 
(we will return to this point). The shift (|19|l leaves the 
measure du invariant, so 



VtK;/3) = I due^i"^-''-"^" X 



(20) 



where the integral has now split into two independent 
factors. The Gaussian integral over du is elementary-42i: 



due-'^^'^"'^-^- 



(2^/3-1 



3A'atom/2 



(det'Z?) 



-1/2 



= Ci/3- 



,/2 



(21) 



where we have used — kT and det' denotes the deter- 
minant without the zero eigenvalues. On the second line 
we note the simple power-law dependence on /3; the con- 
stant factor Ci is ultimately irrelevant. In order to sim- 
plify notation, we define Kjk = {fji^D'^fku) , where 
again a suitable regularization is implied on the right- 
hand side. Then we have 



,/2 



(22) 



Now we compute the next integral, again by introducing 
a suitable shift in the variables: 

Xj = Xj - i(3KjkUk (23) 

so that 

Zpot{uk;l3) = Ci X (24) 
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Again the Gaussian integral is elementary 



where only the dependence on (3 is relevant 44 Finally, we 
have the expression we need: 



-Af„odc)/2 



The CG (harmonic) potential energy is then 



Epotiuk) = -5/3 log 



pot 



(26) 



(27) 



3 1 

= -^{N^tom - Nnodc)kT + -UjKjkUk (28) 

This expression shows that the CG potential energy is 
exactly given by the sum of two contributions. First, 
each degree of freedom that has been integrated out con- 
tributes a thermal energy of ^kT to the potential energy. 
And second, the remaining CG degrees of freedom expe- 
rience a harmonic potential of the form ^ujKjkUk- The 
calculation of the kinetic energy is essentially identical 
in form, but a bit more simple due to the diagonal mass 
matrix. We note that while this derivation is mathe- 
matically elegant, it finesses many subtleties in a way 
that may not give the reader complete confidence in the 
derivation; e.g. we have ignored the fact that the dynam- 
ical matrix is singular and we have not been very careful 
about the values of the irrelevant constants Ci and C2. 
We present a more careful derivation in Appendix IXlthat 
has the distinct benefit that the reason for, and extent 
of, the non-locality of the stiffness matrix is readily ap- 
parent. We now return to the CG energy, and once again 
separate the spatial and nodal indices for the following. 

The full CG energy H16I) for a monatomic harmonic 
solid of A'atom atoms coarse grained to N^odc nodes is 
then given by 

E{Uk,Uk) = Uint + i^jk • Ufc + • KjkUk) , 

(29) 

where the contribution of the internal degrees of freedom 
is 

Uint = iVatomS™'^+3(A^atom-Afnodc)fcT. (30) 

Mjk and Kjk are defined as follows. The mass matrix is 



M. 



jk 



E/^m^^'Am) (31) 
I'^^Nj^Nkn (monatomic), (32) 



where the second line applies to monatomic solids with 
atomic mass m. This may be expressed in matrix nota- 
tion as 

M,k = \{NN^){N Mllj,N^y\NN^)\ (33) 

L J jk 

= m{NN'^)jk (monatomic), (34) 



where the MD mass matrix is M^IP = m^Si^i^. Note that 
each of the quantities in parentheses () in Eq. 133|) is an 
Nnode X Nnode Square matrix. 

At times, it may be desirable to use a diagonal approx- 
imation of the mass matrix, often called the lumped mass 
matrix in the FEM literature. In CGMD it is given hy^ 



M. 



lump 



(35) 



The CGMD stiffness matrix is given formally by 



Kjk 



(E/^a.^a:-'a-) (36) 

= \{NN^){ND^^N^)~^ (NN'^)] , (37) 
L J jk 

where each of the entries on the second line is a matrix. 
The inverses in (I31f) . H33|) . H36|l and (|37|) are matrix in- 
verses. 

Remarkably, there is another way to write the CGMD 
mass and stiffness matrices that does not require two 
inverses. These forms are derived in Appendix ^ and 
the notation is explained there. In particular, they are 
given by 



K,k = N,,D,,Nk. - D^^D-^Dl^ (38) 
M,k = N.^m^Nk^ - M,x AC>4>t (39) 



We do not currently know of an easy way to show the 
equivalence of Eqs. (|36|l and H38() . but we have checked 
that they are numerically equal. Each form has its pre- 
ferred uses. In the calculation of the CGMD spectra 
and any other application in which the stiffness matrix 
in reciprocal space is needed, Eq. is advantageous. 
It is formally simpler, but it suffers from requiring two 
inverses and from formal singularities due to the zero 
modes of the dynamical matrix. Both of these drawbacks 
disappear in reciprocal space, where the dynamical ma- 
trix is diagonal (in the Fourier transform of the nodal 
indices). On the other hand, Eq. H38|) is well defined, 
and it only requires one inverse. Typically, the inverse 
in the atomic indices is taken in reciprocal space making 
use of the perfect crystal space group symmetry, whereas 
the second inverse (|36|l is in the nodal indices for which 
reciprocal space offers an advantage only in special cases 
where the mesh is uniform. For irregular meshes, Eq. 
(|55)l is preferred. 

Consider the form of the second expression for the stiff- 
ness matrix (|38|) . The first term represents a form of 
coarse graining in which each atom is forced to be ex- 
actly at the position defined by the interpolation func- 
tion. Within the context of CGMD, we will refer to the 
approximation where the other terms are neglected as 
the rigid approximation. To be precise. 



Pr^ Rigid Approximation 



(40) 



by which we mean that P^, defined in Eq. ljA7p is set 
to zero in all of the subsequent CGMD formulas. For 
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instance, both Dj^ and Mj^ vanish in the rigid approxi- 
mation, so only the first term survives in Eqs. (|38|l and 

Now let us consider the second term in the stiffness 
matrix H38() . This kind of term has not been discussed in 
the context of concurrent multiscale modeling previously, 
and it is very interesting. It is in the form of a lattice 
Green function contribution to the stiffness. According 
to our principle of microscopic-macroscopic correspon- 
dence, the atomistic degrees of freedom are assumed to be 
a best fit to the coarse-grained degrees of freedom. Even 
at zero temperature, this requirement does not necessar- 
ily mean that the atoms will be position exactly where 
the interpolation function would put them. Instead, they 
typically relax to a lower energy configuration. This re- 
laxation in the short wavelength degrees of freedom in- 
troduces a non-local coupling between the coarse-grained 
degrees of freedom through the Green function . The 
appearance of a Green function in a relaxation problem is 
natural. Consider a system governed by the elastic ener- 
getics E = ^Ufj^D^yUi, — fftUfj,. The minimal energy state 
is Uf_, = D~^f^ with the energy E^^in = -\u^^D^^u^. 
The Green function thus arises naturally in the energy 
of the relaxed state. In CGMD, it is the internal modes 
that can relax, and so it is the Green function of these in- 
ternal modes that enters the CGMD Hamiltonian and in- 
troduces the non-locality. This kind of weak non-locality 
is not present in finite element modeling, but it is en- 
tirely physical. In fact, continuum formulations of non- 
locality in elasticity have been introduced to account for 
size effects in dislocations, crack tips and other nanoscale 
structures It may be possible to neglect the non- 
locality for a particular application, but it is real and it 
arises naturally in CGMD. 

The energy (|29|l contains terms representing the aver- 
age kinetic and potential energies, plus the thermal en- 
ergy term expected from the equipartition theorem for 
the modes that have been integrated out. As mentioned 
above, this Hamiltonian continues to work for polyatomic 
solids, in which the optical modes may be coarse grained 
in various ways to represent different physics. 



B. No Ghost Forces 

One goal of concurrent multiscale modeling is to have 
the atoms in the atomistic region behave as closely as 
possible to the way they would if the atomistics extended 
throughout the system. Deviations from this ideal have 
been termed ghost forces. Some deviation is inevitable, 
but two kinds of deviation have received particular at- 
tention as pathologies. Both are at zero temperature. 
First, if the displacement with respect to the equilibrium 
lattice is zero throughout the system, then the atoms at 
the interface should experience no force. Second, if the 
displacement corresponds to uniform strain and the uni- 
formly strained atomistic system is in equilibrium, then 
the atoms at the interface in the concurrent multiscale 



model should experience no force (See the recent review 
article by Curtin and Miller— ). 

CGMD does not suffer from the first type of ghost 
force, as can be seen immediately from the absence of 
terms linear in Uj in Eq. H29|l . The second type of ghost 
force is also absent from CGMD, provided the strain is 
admissible in the space of shape functions. This property 
should be clear from the construction of CGMD, where 
if the best fit interpolation function reproduces the uni- 
form strain at zero temperature, the delta functions im- 
pose that the CGMD energy agree with the MD energy 
exactly. 

To be more precise, the strained state described by 
the underlying atomistic displacement is admissible if 
there exists a set of nodal displacements such that 

= I]uj^j(2^mo) (41) 

j 

= u^fj.Nj^ (42) 

= P^f^. (43) 

where we have used Eq. H1U|) in the second line to express 
the best-fit u., in terms of u^,. The matrix P!^,P , defined 
in Appendix lAi in Eq. (|A6|) . must act as the identity on 
Ui^. We assume that the uniformly strained atomistic 
system is in equilibrium, and hence D^^Ui, — 0. We cal- 
culate the force on node i (which may be an atom at the 
interface) using Eqs. (|38|l and H1U|) as 

- K,j^3 = ~N.,^D^,N,, (/,,u,) - DlD-^^D^ {f,,m 

= N,^,D^,vi, - D^^b-lP^^DruVi, (46) 
= (47) 

Thus the atoms (and nodes) of the coarse-grained sys- 
tem experience no force. In going from the second to 
the third line we have used the admissibility condition 
that P^^ act as the identity matrix on Up, and in the 
following line we used that is in equilibrium and so 
Dt^Mu = 0. This derivation proves that any admissible 
equilibrium atomistic configuration is also an equilibrium 
CGMD configuration. The derivation continues to work 
once anharmonic forces are included, as discussed below, 
since the terms in the CGMD energy up to second order 
in displacements are the same. Thus, CGMD is free of 
ghost forces in both senses of the term. 

C. Anharmonic Lattices 

We have formulated CGMD for an underlying an- 
harmonic Hamiltonian in perturbation theory, assuming 
again negligible diffusion in the CG region. With an an- 
harmonic potential the higher frequency modes compris- 
ing the heat bath do not decouple, and they introduce 
temperature-dependent effects such as thermal expansion 
and thermal softening of the elastic constants. The basic 
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idea of how anharmonicity is treated was presented in 
Ref. |3- The details are presented here for the first time. 

In this subsection we develop a formal analysis of the 
contributions of the anharmonic interatomic forces to the 
CGMD energy and equations of motion. This analysis 
provides insight into how effects in the coarse-grained 
system are linked to their atomistic origins through ana- 
lytic formulas, and hence are very powerful. On the other 
hand, the formulas are sufficiently complicated that a dif- 
ferent approach is employed in practice, and we stress 
this point. The formal developments that follow take the 
perfect T = OK lattice as the reference state; in practice, 
it is much more useful to take the lattice at the temper- 
ature of interest to be the reference state. In that case, 
the harmonic theory presented above may be used, with 
anharmonic effects from the lattice entering into the ref- 
erence state in a quasi-harmonic approach. The thermal 
expansion of the lattice and the thermal softening of the 
dynamical matrix capture the anharmonic forces The 
finite temperature lattice constant and dynamical ma- 
trix are inputs to the CGMD formulation, precomputed 
in conventional MD calculations. While that approach 
is very effective in practice, it is more satisfying from a 
theoretical point of view to have a direct, analytic the- 
ory of the thermal effects in CGMD, and it is to this 
development we turn now. 

The CGMD energy is computed by perturbation 
theory about the harmonic Hamiltonian H15|l . Specifi- 
cally, 



Then 



E = 
Hh = 

H' = 



Ut,t - dfj log J du^dp^ e'"'^ e-"' A, (53) 
E||+E^"M-^..ii., (54) 

/J. ^.V 
OQ ^ 

E E ■ ■ ■ (55) 



n— 3 



H 



MD 



Hh + H' 



(48) 



where tz/j^j is the internal energy for the CG harmonic 
lattice (123 ■ The temperature dependence is now in C/j^^, 
A, and in the perturbation H', but not in Hh- 

To evaluate this integral, the exponential involving H' 
is expanded in a Taylor series in u and/or kT. The re- 
sulting integrals are elementary. It is common practice 
to use a diagrammatic representation of the integrals to 
facilitate bookkeeping in this kind of expansion^ The 
log in Eq. I|53() is then produced by restricting to con- 
nected graphs.^* Similar perturbative approaches have 
been used in many contexts; the application most rele- 
vant to the current study is the Self-Consistent Phonon 
Approximation used in lattice dynamics42iS 

The resulting form of the CG energy for the anhar- 
monic lattice is 

j-k 

CO ^ 

E;;TE^::;r(^)";^•■";:(56) 



where Hh is the harmonic Hamiltonian (|15|l and H' con- 
sists of the anharmonic corrections, which are assumed 
to be small. This is a good approximation in silicon, 
and many other materials of interest below their melt- 
ing point (i.e. at low homologous temperatures). The 
perturbation may be written explicitly as 



H' 



E 



(49) 



where the /x's label the atoms, as before, and the a's label 
components of the vectors. The higher order dynamical 
matrices D^^ 'I^Z taken to be completely symmet- 

ric in their indices from the form of Eq. H49|l . Below we 
occasionally use the notation to represent -DJJj '^" 
schematically in places where the additional concision 
should not lead to confusion. This perturbation theory 
is a low temperature expansion, as seen by switching to 
the scaled coordinates 



A, 



VkT' 



A' 



/kT' 



(50) 

(51) 
(52) 



where now C/jnt is a complicated function of f3, as are 
the stiffness coefficients. Our goal is to calculate this ex- 
pansion analytically at each order of perturbation theory. 
Since the diagrammatic approach may not be familiar to 
all readers, we will not use it for the derivations, but 
we return to it below. For now, consider the integral 
needed to calculate the CGMD energy up to second or- 
der in D^"^^ and first order in D'^^\ These first few terms 
in the perturbation series will capture the leading ther- 
mal and non-linear effects. Higher order terms could be 
calculated similarly, if needed. To this order, the CGMD 
potential energy expansion is 



_!_ r^aia2a3 ~ai -a^ -03 

1 ( r)aia2Q3 „~,ai -02 .-.as 

2 A'1M2P3"a'i"a'2"m3 



(57) 



a-1 r-)aia2a3a4 iT^i --02 -03 -04 . 

I'' -^MlP2M3A'4 "pi "^2 "^3 "/i4 



-Ho 



where the polynomial terms have resulted from expand- 
ing the anharmonic exponential (|55|) in a Taylor series. 
In classic lattice dynamics analysis, a similar expansion 
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is used. There the cubic term is odd in u and integrates 
to zero, while the other two terms give non-zero contribu- 
tions to the heat capacity of the order kT^ In CGMD, 
each of the three terms gives a non-zero contribution, as 
we will now show. 

Integrals of the form H57|l , a polynomial times a Gaus- 
sian, may be evaluated in closed form. The algebra can 
be tedious, so often generating functions are used to sim- 
plify the bookkeeping. The CGMD generating function 
is defined to be 

Zh[J] = Zki„ydii^e-5U.. A..a. e^M-u.A (58) 

where Z^i^ is the part of the partition function coming 
from the kinetic energy and including the nodal velocity 
constraints. The subscript on Zh indicates that it is the 
generating function corresponding to the harmonic the- 
ory. Generating functions are useful because derivatives 
with respect to bring down factors of u^, 

^Ja^" ['^]|,7=0 = ^kin j du^e-^"--^-^"'^ UA A (59) 

so they provide a simple way of calculating the integrals 
of Gaussians multiplied by a polynomial prefactor. The 
CGMD internal energy 1)5 7|l may therefore be computed 
from the generating function as 

U = Ut^, - dp log { exp [~f3H' {dj)] Zn [J] \ } (60) 

where the factor exp[—f3H'{dj)] is intended to be ex- 
panded in its Taylor series for practical calculations. 

We return to the anharmonic terms below [cf. Eq. H71|l ]. 
but first we compute the partition function. In order to 
compute Zh [J] we can complete the square and relate the 
result to Zh[0]: 

Zh[J] = Zki„e5J.-^M.3J. y"dW^e^5w.-i3..w„ ^ (g^) 

^ ^h~^.-D-li. z,,[J = 0; u, - Jj^D-l ■ J,](62) 
= Z,JJ = 0]e-3JM-GM--J.^"H,-J« (63) 

where = — ■ 3^. We have used Eq. (^5)1 to 
derive the third line. The Green function G^^ and the 
external field are given by 

G,. = D-'j,,K,t,hxDll - D-l (64) 
= u,K,kfk.D-^ (65) 

with = /3i/2H^. 

Using the projection matrices defined in Appendix 1X1 
we can rewrite this expression in a more transparent 
form. The projection matrices are given by = 
Nj^ifju and P^^ = - P'^^ . Using them together with 
the formula for Kij (|36|l we have the alternate formulas 

G,. - P^p {D-p^finK.kfkxDll - D-^) P.i (66) 
= u,N,^ + xx,K,kfk.D-lPt^. (67) 



where D^^ is defined by Eq. (|A16|l i^ These formulas 
show that in the rigid approximation G^i, — and 

= UjNj^. The result is that each occurrence of 
in the MD potential energy is just replaced by Nj^Uj 
in R-CGMD. Relaxation of the unconstrained degrees of 
freedom makes the further contributions involving P^. 

To be more explicit, the rigid approximation to CGMD 
freezes the unconstrained (unresolved) modes. Formally, 
this is accomplished by setting the orthogonal projectors 
to zero: 0. In this approximation the generating 

function simplifies 

ZH[J]U^,^^Zn[Q\e-'P'''^^^^^^^ (68) 

where all of the Green function contributions have been 
eliminated. Similarly, the mass and stiffness matrices in 
Z/i[0] involve only the first terms in Eqs. (|38l) and l|39|l . 
The Rigid CGMD (R-CGMD) Hamihonian is then given 
by 

£;(ufe, life) = J/int + u, • lij -I- UMD{Nj^,\ij), (69) 

where Umd is the full, non- linear MD potential energy 
evaluated with = Nj^^Uj, and My = Nj^rrifj^Nkfj.- The 
rigid approximation neglects thermal effects and relax- 
ation of the internal degrees of freedom (i.e., those that 
have been integrated out). This expression is a conserved 
energy, and it is free from spurious forces commonly re- 
ferred to as 'ghost forces' in discussions of concurrent 
multiscale modelingiiSS It is also as computationally ex- 
pensive as MD, because unlike the stiffness matrix, the 
non-linear potential cannot be pre-computed except per- 
haps for some toy potentials. Additional approximations 
are needed. One set of approximations based on the use 
of representative atoms and the Cauchy-Born rule leads 
to the Quasicontinuum Methodi^ Thus, the derivation 
above shows the relationship between CGMD and the 
Quasicontinuum Method, in particular that the Quasi- 
continuum Method is closely related to the zero temper- 
ature rigid approximation of CGMD. Another success- 
ful concurrent multiscale technique. Coupling of Length 
Scale (CLS)^^, can also be related to CGMD. It involves 
taking the large- iV asymptotic limit of the stiffness ma- 
trix in the rigid approximation, where N in this case is 
the number of atoms per element. The large- iV limit 
yields continuum mechanics within the element; i.e., a 
finite element representation of continuum mechanics. 
Both CLS and the Quasicontinuum Method make ad- 
ditional approximations at the MD/CG interface. 

We now consider the physics of the terms in the CGMD 
energy beyond the rigid approximation. These are ef- 
fects that arise due to the anharmonicity of the inter- 
atomic potential, and include both thermal effects and 
non-linear relaxation. The diagrammatic approach is a 
convenient and powerful way to analyze perturbation the- 
ory at higher order, such as the anharmonic contributions 
to CGMD. The quantities represented by the diagrams 
are typically generated from a relatively simple set of 
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FIG. 2: Feynman diagrams of connected graphs contribut- 
ing to the anharmonic CGMD internal energy. The thick 
external legs represent factors of H, which in the usual termi- 
nology correspond to dressed external fields. Physically, they 
are contributions from the CG fields with some zero tempera- 
ture relaxation of the internal degrees of freedom. The dashed 
internal lines represent factors of G that arise from the fluctu- 
ations of the internal degrees of freedom. These fluctuations 
are temperature dependent, and each factor of G comes with 
a factor of kT). 



denote contracted tensor indices. 

At a given level of perturbation theory, the graphs pic- 
torially represent the hierarchy of contributions to the 
CGMD energy. Although graphs that split into more 
than one disconnected sub-graph contribute to the par- 
tition function, only connected graphs contribute to the 
free energyii^ These are shown in Fig. [3 Consider first 
the graphs with no external legs. These terms in the 
perturbation theory are independent of u^; they only 
contribute to the internal energy resulting from degrees 
of freedom that have been integrated out. Thus they 
build up the non-trivial temperature dependence of Umt- 
The graphs at the other extreme-in particular those that 
have no dashed lines and are therefore independent of 
G^i/-make up the rigid approximation that was discussed 
above, but now including some internal relaxation. 

Using the generating function (|63|l . we have calculated 
the CGMD internal energy (|57() up to first order in D^^^ 
and second order in D^^^: 



Z = 



-1/2 



T~jaia2a3 Q „ a _ _ 



(70) 
(71) 



tl mi 

72 

■ -'^Mi ^2^3 M4 c'j^; "jii ^Jii 



^TuTuXZlldj^i dj^3 dj^s dja, aj»5 dj^ 



Zh[J] 



.7=0 



where Zh[J] is given by Eq. (|63|l . The derivatives act on 
the Gaussian generating function 163(1 to give 



rules, known as Feynman rules. Due to space limitations 
we cannot provide a thorough review of the rich math- 
ematical structure encoded in Feynman diagrams, but 
direct the interested reader to one of the numerous texts 
on the subjectf^ We have derived the Feynman rules 
for CGMD, but they will not be presented here since 
they are not actually needed. Instead, we will employ 
the set of Feynman diagrams as a pedagogical device to 
consider the form of various contributions to the CGMD 
energy. The Feynman diagrams up to second order in 
D'^'^^ and first order in D'^^^ are shown in Fig. |21 In the 
diagrams, the thick external legs represent factors of H^; 
i.e., factors of the CG fields Uj, suitably dressed to ac- 
count for elastic relaxation of the internal degrees of free- 
dom. These factors are playing the role of external fields 
in the statistical field theory of the equilibrium state of 
the internal degrees of freedom. The thinner, dashed 
lines represent factors of the Green function G^^. The 
Green function accounts for thermal fluctuations of the 
internal degrees of freedom, and indeed each factor of 
G^i/ brings with it a factor of kT. Vertices in the graphs 
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Next we find the Helmholtz free energy F = —kTlogZ: 



F = Fh. 



armonic 



H'{n^ = H^) 
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(73) 



where ^harmonic IS the CGMD free energy for a harmonic 
crystal [cf. Eq. 



armonic 



The free energy is suitable for isothermal CGMD simu- 
lations. On the other hand, the internal energy is appro- 
priate for adiabatic simulations, and hence more closely 
related to MD simulations without a thermostat. It is 
given hy U — dp {(iF): 



U ^ 



-^harmonic H" ^ ('-^M " -^m) 

1 
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(kT\^ rioi0203 T^aiaaa^ ^0104 ^a^a^ ^a^ag 
V / M1M2M3 M4M5Me M1M4 M2MS M3M6 



(kT)^ r)Oi02a3 jjaiasae 7^0104/^0203/^0508 
g V / M1M2M3 M4MSM6 M1M4 M2M3 MsMe 

+ . . . (76) 

where ^harmonic is given by Eq. (|29|l . Now T is under- 
stood to be a function of the (constant) entropy and 
the state of deformation, although for many applications 
since the deformation in the CG region is small, it can 
be treated as constant. For reference, the expressions for 
the Green function G^^?,^ and the field H"; are given in 

M1M2 M 

Eqs. (|dd)1 and l|67|) . respectively. 

We could continue to calculate the terms in the GGMD 
Hamiltonian to higher order in the MD anharmonic cor- 
rections and/or the CGMD thermal perturbation expan- 
sion. The number of terms grows rapidly, and we quickly 
reach the point of marginal returns; i.e., the point at 



which the added complexity is no longer rewarded with a 
commensurate improvement in accuracy. However, there 
are certain kinds of contributions where improvements 
are possible. In particular, it is important to capture the 
first non-trivial effects in the expansion. The terms cal- 
culated thus far have contributed to the harmonic Hamil- 
tonian, the zero-temperature anharmonic terms and the 
energy of the internal modes. We have also calculated the 
leading contributions in the free energy to the thermal 
expansion and temperature-dependence of the stiffness, 
i.e., the terms proportional to and Uj^ (g) Uj^. In the 
internal energy, these contributions are higher order in 
the anharmonic lattice expansion, so we have not calcu- 
lated them yet. Of course, there are thermal corrections 
to the higher order stiffnesses, as well, but they are not 
as important. 

We conclude this section with the calculation of the 
leading temperature-dependent, quasi-harmonic contri- 
butions to the CGMD Hamiltonian (internal energy). 
The graphs contributing the leading temperature depen- 
dence to the one-point function (governing thermal ex- 
pansion) are shown in Fig. 13 The term "one-point" 
means a single leg, and hence a single factor of Uj. Con- 
sider the lowest order term in kT in Fig. |31 as it enters 
the free energy: 
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where we define the temperature-dependent fj' (not to be 
confused with /j-^) to be: 



- 2 ^MiM2M3 G;\Z h.D-l^K,k (78) 



where we have used the expression (|65f) for . We have 
written the leading one-point term from the free energy 
in Eq. (|77(l in a form familiar from finite element models 
as an eigenstrain (cf. Section 2.12 of Hughes^i). As in 
conventional finite element formulations in uniform tem- 
perature, this term is a total difference (the discrete ana- 
log of a total derivative), and it only enters through the 
boundaries of the simulation. In a homogeneous system 
with free surfaces, the result is expansion of the system 
as the temperature increases. The expansion is propor- 
tional to kT and D^'^^ to leading order, as expected based 
on conventional lattice dynamicsi^ Note that since the 
linear term is a total difference, the discussion of ghost 
forces given above continues to be valid in the case of 
thermal expansion. 

The graphs contributing the leading temperature de- 
pendence to the two-point function governing the tem- 
perature dependence of the adiabatic elastic stiffness are 
shown in Fig.^ The leading contributions to the isother- 
mal stiffness are the three graphs containing Gi but with 
Gi replaced by G. As expected, the first thermal contri- 
butions to the adiabatic stiffness are of order (fcT)^ due to 
the third law of thermodynamics (the Nernst theorem), 
whereas the first thermal contributions to the isothermal 
stiffness are of order kT. 
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FIG. 3: The leading anharmonic contributions to the parti- 
tion function afTecting the thermal expansion in CGMD. 



Finally, we reemphasize that we have calculated prop- 
erties within the CGMD perturbation theory to show 
that the theory is consistent and to gain some theoreti- 
cal understanding of the interplay of anharmonicity and 
temperature in coarse-grained systems. In practice, we 
take the reference configuration to be the crystal lattice 
at a particular temperature with the corresponding finite 
temperature dynamical matrix. It would be a tautology 
to compute thermal expansion or thermal softening in 
CGMD: they agree with the MD result identically by 
construction. 
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FIG. 4: The anharmonic contributions to the partition func- 
tion affecting the temperature dependence of the adiabatic 
CGMD stiffness. The Green function Gi is represented by a 
dashed line with a large filled circle. It is defined in the upper 
box. 



V. THERMAL FLUCTUATIONS 

In the previous Section we developed a description of 
the average motion of the collective degrees of freedom. 
We have shown that the short wavelength modes con- 
tribute to the mean motion, even though they are not 
explicitly present. They have the additional effect of 
inducing small fluctuations about this mean. The CG 
modes behave as if they were in a Brownian heat bath; 
they are jostled by small, random interactions with the 
(invisible, missing) short wavelength modes. In this sec- 
tion, we analyze these interactions in detail. 

Consider the fluctuations of the CG fields in harmonic 
MD in thermal and mechanical equilibrium, 

(ujUfc) /,^/fe,Z-iy"dx^dp^ 6-'^^"° u^u, (79) 

= fjJkukTdD.JrlogiD) (80) 
= kTf,^D-lhu (81) 
= kTKJ^' (82) 

This calculation shows that the two point isothermal cor- 
relation function of the CG displacements grows linearly 
with the temperature. The result is reasonable and could 
equally well be understood in terms of the equipartition 
theorem. As the temperature increases, the average po- 
tential energy increases and so the mean amplitude of the 



harmonic oscillations increases. Note that the amplitude 
decreases as the stiffness Kjk increases. 

We can repeat the calculation in harmonic CGMD 
again in thermal equilibrium for comparison: 

(ujUfc)T ^ Z-^ Jdujduj e-'^"'^'^'"'' UjUk (83) 
= fc^if-i (84) 

where Hcgmd is given by Eq. I|29() . The result is exactly 
the same as the correlation function of the CG fields in 
MD. This equivalence of the correlation functions holds 
for higher correlation functions as well, due to Wick's 
theorem.^- 

While the equilibrium properties agree, there are dif- 
ferences in the time autocorrelation functions of the dis- 
placement. These differences result from the influence of 
the short wavelength modes that are not represented on 
the mesh. These additional modes acts as a heat bath, 
exerting random and dissipative forces on the CG fields. 
The effect is entirely analogous to Brownian motion, in 
which a large particle is jostled by the thermal motion of 
the unseen atoms in a liquid surrounding itS 

There are several practical issues that arise in the anal- 
ysis of CGMD simulations regarding fluctuations. In con- 
ventional MD, the mean atomic kinetic energy is used to 
compute the temperature. We now investigate whether 
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a similar connection holds in the CG region of CGMD 
(naturally it continues to hold in the MD region if it is 
in thermal equilibrium). We consider the simplified case 
in which there is no large-scale motion such as center-of- 
mass motion; our derivation continues to hold provided 
any such non-equilibrated modes are subtracted from the 
nodal velocities prior to analysis. Consider the mean- 
squared velocity of node i in thermal equilibrium. We 
use the partition function for CGMD to calculate the 
expectation value of \ui\^ in the canonical ensemble at 
thermal equilibrium: 



= Z, 



kin 



dii 



J=0 



j=o 



(85) 
(86) 
(87) 



the form of random and dissipative forces. The form of 
these generalized Langevin forces may be computed using 
statistical mechanical (Zwanzig-Mori) projection opera- 
tor techniques, although it is beyond the scope of this 
article. The resulting spatio-temporal memory kernel 
has been described elsewhere^ The random, dissipative 
forces not only act to bring the CG degrees of freedom 
into equilibrium with the internal degrees of freedom, but 
they act to absorb short wavelength modes incident on 
an interface where the mesh is refined. In principle, they 
also include the propagators that reconstruct waves on 
the far side of the CG region if the mesh is refined again; 
in practice, these propagators appear to be very expen- 
sive to implement computationally. In fact, one of the 
challenges of memory kernels is their computational ex- 
pense both in terms of the memory required to store the 
recent history and in terms of the demands they place on 
parallelization to make the code suitable for supercom- 
puters. 



where the factor of 3 in the final line is the number of di- 
mensions. The potential energy contribution to the par- 
tition function is irrelevant to this calculation and has 
been suppressed. The important result expressed in Eq. 

is that the mean-squared velocity is directly propor- 
tional to the temperature, just as it is in conventional 
MD. 

In fact, the MD result is recovered by replacing Ma 
with the mass of atom i. This calculation also applies 
to concurrent multiscale models that use a lumped mass 
matrix. For CGMD with a non-diagonal mass matrix, it 
is only the diagonal of the mass matrix that affects the 
amplitude of oscillation of a given node. The off-diagonal 
terms do introduce correlations between the velocities of 
neighboring nodes that would not be present for a di- 
agonal mass matrix. Because of these correlations, the 
simple atomistic relationship between the mean-squared 
velocity and the temperature that follows directly from 
the equipartition theorem must be modified according to 
Eq. ^ for use in CGMD. 

Thus far in this article we have treated CGMD as a 
system that conserves energy, and these random, dissipa- 
tive forces are absent. In particular, the evolution of the 
system generated by the Hamiltonian ||2J) conserves the 
energy given by that Hamiltonian. Thermostats might 
be added to simulate the electron-phonon coupling; i.e., 
the interaction of the lattice vibrations with the elec- 
tronic degrees of freedom. Such additions violate energy 
conservation?^ since energy can flow to and from the 
heat bath and the system becomes an open system. 

Even neglecting the electron-phonon coupling, the 
coarse-grained system of solid mechanics described by 
CGMD is an open system. In the full MD system en- 
ergy can flow between modes that would be retained in 
CGMD and those that would be integrated outi^ This 
interaction implies that the CGMD energy is conserved 
only on average, and that additional interactions are 
present in reality. These additional interactions take 



VI. IMPLEMENTATION DETAILS 

In practice, CGMD is run much the way conventional 
MD would be. The forces in the CG region are deter- 
mined by the CGMD stiffness matrix and the nodal dis- 
placements; the forces in the MD region could be deter- 
mined this way, as well, but since we have shown that 
the forces in the MD region are just the usual MD forces, 
the full MD potential is used to calculate the MD forces. 
Using the accelerations, the velocity Verlet time integra- 
tor is used to evolve the system in time»^ The same 
time step is used throughout the simulation. In princi- 
ple, the natural frequency in the CG region is lower as 
the mesh size increases, and a longer time step could be 
used there; in practice, the CG region entails relatively 
little computational expense, and there is little motiva- 
tion to introduce a spatially varying time step that could 
cause subtle problems. 

One difference from MD and conventional FEM is that 
the topology of the CG mesh is not allowed to change. 
Neighboring nodes remain neighboring nodes throughout 
the simulation. The topology of the mesh is determined 
by a cell list, which contains the nodes associated with 
each cell in the mesh, and a face list, which lists pairs 
of neighboring cells. An edge list, which lists pairs of 
neighboring nodes, is also generated. 

The stiffness matrix Kij is to be computed once at the 
start of a simulation, and it remains unaltered during the 
subsequent dynamics. It does not matter whether atoms 
vibrate across cell boundaries, as long as the crystal lat- 
tice topology does not change and diffusion is negligible. 
The elements of Kij decrease exponentially with distance 
from the diagonal, and in practice it is necessary to trun- 
cate the stiffness matrix in order to control the memory 
and CPU requirements for simulating large systems with 
irregular meshes. Both requirements scales as at least 
0{N'^) if the full stiffness matrix is retained, and this 
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scaling can be reduced to 0{N) if matrix elements are 
discarded if Kij < eKn for some small number e. This 
approach allows the simulation of billion atom systems 
(greatly coarse-grained) on desktop workstations with- 
out approximation beyond those presented here. 

For T ^ the finite temperature dynamical matrix 
should be used for U^j/. This quasi-harmonic approx- 
imation ensures a consistent thermodynamics, and it 
effectively sums the two-legged diagrams of the finite- 
temperature anharmonic perturbation theory (i.e., those 
terms second order in the CG displacement). For exam- 
ple, in ergodic systems the time average of the kinetic 
energy term in the CG energy is related to the tem- 
perature through a equipartition expression. In general, 
the dynamical matrix may depend on other macroscopic 
parameters, as well, such as slowly varying external mag- 
netic and electric fields. should be evaluated under 
these conditions. Also note that while the harmonic ap- 
proximation may be good in peripheral regions, it may 
not be appropriate for the important regions. We have 
shown that the CGMD and MD equations of motion 
agree in regions where the mesh coincides with the atomic 
sites. In these regions, the full MD potential may, and 
should, be employed, so that effects such as diffusion and 
dislocation dynamics are allowed. 



A. Normal Modes and the Inverse of D^i, 

In order to simulate large CG regions, it is necessary to 
take some measures to increase computational efficiency. 
One such trick is to make use of the long-range order in a 
single crystal to facilitate the computation of the stiffness 
matrix. The eigenstates of the dynamical matrix D^^ are 
plane waves. For monatomic lattices they correspond to 
the normal modes of the system, the longitudinal and 
transverse phonons in the acoustic and optical branches 
that are familiar from lattice dynamicsi^SiS In reciprocal 
space, where the basis elements are exactly these plane 
waves, the dynamical matrix is diagonal. The inverse 
of the dynamical matrix is then trivial to compute, and 
the subtlety of inverting a singular matrix is eliminated 
because reciprocal space naturally factorizes into a direct 
product of the three zero modes with k = and all of 
the other modes with non-zero eigenvalues. 

To be specific, the equation for the stiffness matrix H37|) 
becomes 

= {NN^),^[N^{--k)D-,\-k)N,,{-k)]'\NN^m 

where k = is explicitly omitted from the sum. Now the 
inner matrix inverse is the inverse of a 3x3 matrix. 



B. Shape Functions 

This section discusses particular choices for interpo- 
lation functions. Compatible combinations of these are 



also allowed, as in FEM with wedge and brick elements, 
for example. We emphasize how different choices for in- 
terpolation functions meet the requirement of meshing 
the crystal lattice at the MD/CG interface. The usual 
linear interpolation functions for tetrahedral elements^- 
are the simplest functions meeting the three criteria (|IJ 
- (Im)) in Section UTTl They are defined such that Nj{x.) 
is 1 at node x^, it goes linearly to zero at the nearest- 
neighbor nodes, and it vanishes outside of the nearest 
cells. Suppose x is in the kth element with nodes x/j. 
where j — 1, . . . , 4. Then the interpolation functions are 
given by the volumnal or natural tetrahedral coordinates: 



X • logyfc(xfei, . . . ,XfeJ 
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(90) 
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where = {xk,yk, Zk) and we have written the tetra- 
hedral volume as a determinant. The interpolation func- 
tion for node kj is simply the volume of the tetrahedron 
formed by x and the other three nodes divided by the 
volume of the entire tetrahedral cell. These functions 
are clearly C° continuous and independent. It is eas- 
ily checked that they are linear and form a partition of 
unity. They also have the desirable properties of locality 
and ease of use. The locality property is particularly im- 
portant for our applications, since the domains requiring 
an atomistic treatment are localized to small regions of 
the system. 

Another basis we have found useful is the set of the 
longest wavelength normal modes. These functions sat- 
isfy the less stringent basis properties: (i') linear indepen- 
dence, detA'j(xfc) ^ 0, and (ii') representation of unity, 
1 = ^2^=1^' Nj (x) for some constants cj . This basis 
provides a check of the CG Hamiltonian (j2Jl, since these 
functions are the optimal basis for a regular CG mesh — 
the phonon spectrum comes out exactly correct, apart 
from the missing short-wavelength modes. The disad- 
vantage of this basis for irregular meshes is that it is 
nonlocal and the short wavelength modes that should be 
supported on the finer parts of the mesh are absent. In 
particular, the stiffness matrix elements decrease as 



1 



(92) 



instead of decreasing exponentially with distance in the 
local basis case. The short wavelength modes could be 
restored locally through the use of a wavelet basis, in 
principle, but we have not implemented a wavelet-based 
version of CGMD. 

In many cases, higher order polynomial interpolating 
functions are the basis of choice. Generalizations of the 
8-node brick used for hexahedral lattices^ are particu- 
larly easy to implement. For example, a generalized 8- 
node brick is the element we used to calculate the CGMD 
spectra for solid argon and tantalum presented in Section 
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IVII Al below. First, consider a simple cubic lattice. The 
basic 8-node brick involves interpolations functions of the 
form^ 

A^r'^'lCa) - i(l + a) (1 + 6) (1 + 6) (93) 

where throughout this section we use the scaled coordi- 
nates G [^1; !)• ^1 is associated with the corner node 
at ^ = (1,1,1). The other 7 interpolation functions are 
generated by the action of the point group Oh on Ni . 

The cubic elements may be applied to a variety of 
hexahedral elements by mapping the real space coordi- 
nates onto the [—1, 1)'^ cube in the coordinate space us- 
ing the standard multi-linear coordinate transformation 
often called "natural coordinates" for the hexahedron in 
the finite element literature^ For our purposes, some es- 
pecially important cases are the monatomic Bravais lat- 
tices, such as face-centered cubic (fee) and body-centered 
cubic (bcc) lattices. Suppose are the basis vectors in 
real space, and are the reciprocal basis vectors such 
that a^j • bfe = Sab- For example, in the fee lattice the 
basis vectors could be chosen to be ai = (0, ^, ^)a, a2 = 
(^,0,5)0 and a3 = (i,i,0)a, and bi = (— l,l,l)/a, 
b2 = (1, —1, 1) /a and ba = (1, 1, — l)/a, where a is the 
lattice constant. Then interpolation functions on the 
Bravais lattice are given by 

^Bravais = N]'''°'%£_a = 2ba • X - 1) (94) 

using the shape functions defined in Eq. 1)93(1 . Shape 
functions for the bcc lattice can be constructed in the 
same way, with the basis vectors be chosen to be ai = 
(i,i,-5)a, a.2 = (5,-^,^)0 and ag = (-i,i,i)a, and 
bi = (l,l,0)/a, b2 = (l,0,l)/a and bg = (0,l,l)/a, 
where a is the lattice constant. 

One drawback of the fee and bcc shape functions H94() 
is that they break the point group symmetry of the lat- 
tice. Acting on the mesh with an element of the point 
group returns a new mesh with the same nodes but often 
a different set of cell boundaries. One example is C|, the 
90° rotation about z. It changes the cell boundaries, as 
can be seen by its action on i(ai -f a2) [i.e., ^ = (i, ^,0)]. 

C| maps this face point to ^' = (1, 0, — |), a point on an 
edge of the original mesh. Another way to understand 
this symmetry breaking is that we made a choice when 
we selected the basis ai = (0, i, i)a, a2 = (i,0, i)a 
and aa = (i, i,0)a. Had we selected another basis, say 
a'l = (-i,0,i)a, a'2 = (0, i,i)aand ai^ = (-1,1, 0)a, 
then the mesh would have been different. It would have 
different cell edges and faces, even though the cell nodes 
would be the same. More importantly for our purposes, 
a displacement field interpolated using one set of shape 
functions cannot, apart from a few special cases, be rep- 
resented exactly with the rotated shape functions. The 
results of CGMD modeling then depend to some extent 
on the choice of basis and associated shape functions. 

The symmetry breaking is small and of no consequence 
in most applications; however, we are interested in us- 



ing wave spectra as a test of CGMD, plotting the spec- 
tra along high symmetry directions in the Brillouin zone. 
The symmetry breaking is somewhat troublesome in this 
case because high symmetry directions no longer possess 
the high symmetry and directions that are supposed to 
be equivalent by symmetry are not. We have developed 
a symmetrization procedure to eliminate the effects com- 
pletely. Its use is limited to applications where high sym- 
metry is important, such as wave spectra, so we present 
it below where the spectra are calculated. 

Another approach to this problem is to introduce poly- 
nomial bases that respect the point group symmetry. The 
well-known serendipity functions^iiS are an example of a 
minimal polynomial basis that respects cubic symmetry. 
The serendipity functions may be generalized in a way 
that makes them suitable for a cubic fee or bcc cell such 
that there is a node for each atom in the cubic unit cell, 
so that the MD degrees of freedom are recovered in the 
atomic limit. These are different than the usual serendip- 
ity functions and, indeed, are not as well suited for con- 
ventional FEM applications because of the location of the 
nodes (e.g. the bcc element has an internal node)i^ They 
do, however, meet our need to match the atomic lattice 
and preserve symmetries in the atomic limit. For exam- 
ple, the face-centered cubic (fee) serendipity functions 
have nodes at the corners and the at the middle of the 
faces of a cube. To the best of our knowledge, this kind 
of FEM interpolation function has not been used previ- 
ously. They are well-suited to Bravais lattices, because 
of their simple action under the point group symmetry. 
Consider a cubic unit cell of an fee lattice with local co- 
ordinates £,a, where — 1 < < 1 for a = 1,2,3. Define 
the function 

iVl(ea) = ^(l+Cl)(l+6)(l+e3) X (95) 

[2(^1 + 6 + 6 - 1) - (66 + 66 + 66)] 

associated with the corner node at ^ = (1,1,1) and the 
function 

A^9(6) = ^(i + 6)(i-6')(i-6') (96) 

associated with the face node at ^ = (1,0,0). Basis func- 
tions, Nj{£^a) associated with the other nodes are gener- 
ated from H95|) and (|96|l by the action of the point group. 
These functions satisfy the strong requirements for an in- 
terpolation basis. Each function vanishes outside of the 
cells containing the corresponding node, and it goes to 
zero at the opposite faces of those cells: it is local and 
continuous. Also, taken together they form a partition 
of unity, as is easily checked. These 14 functions com- 
prise only part of the set of 26 polynomials of order at 
most (2,2,1); i.e. the set of polynomials with terms no 
higher than x^y^z (or x^yz"^, etc.). But they are speci- 
fied uniquely by the three basis requirements and the fact 
that they respect the point group; i.e. a point group oper- 
ation which leaves a particular node invariant also leaves 
the corresponding function invariant. These fee shape 
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functions are most useful for testing purposes such as the 
computation of phonon spectra and scattering properties 
where it is desirable to maintain as many symmetries as 
possible. 

Other (novel) fee symmetric bases are possible if one 
abandons the notion of self-contained elements. Typi- 
cally, equations of motion in finite elements are assembled 
element by element. In CGMD, the interaction between 
two nodes decreases exponentially with their separation. 
Since the interactions are not contained within an ele- 
ment (in fact there is no absolute cutoff to their range in 
principle), the equations of motion are not constructed 
element by element, and the role of the elements is simply 
to guide in the construction of an interpolation basis. So 
we can consider an fee lattice of nodes as four interlaced 
simple cubic lattices. The interpolation function for the 
corner node at ^ = (1,1,1) is given by 

Nii^a) - j^(i + a) (1 + 6) (1 + 6) (a + 6+6-1) 

(97) 

and the functions for the other corners follow from sym- 
metry. This completely determines a basis set which 
satisfies the criteria of locality and continuity. It does 
not satisfy the partition of unity requirement in the 
strictest sense, since the uniform displacement mode is 
over-represented: it is represented once for each sublat- 
tice. A constraint must be introduced that the uniform 
displacement on each sublattice is equal to the mean dis- 
placement: 



u = ^ Uj /Nnode for all sublattices p. 



(98) 



Note that all the nodes are equivalent. This is possible 
when the nodes are in a Bravais lattice such as fee, but 
it is not true of the fee basis in H95() and H96|l . 

Interpolation functions for other crystal lattices are 
also available, either because they exist already in the 
finite element literature or because they are easily gener- 
ated. We consider a few cases here. 

Interpolation functions for the bcc lattice may be con- 
structed in a similar fashion. They are of the order 
(2,2,2). In particular, the shape function for the center 
node (0,0,0) is 



^9(6) = (1 - 6') (1 - el) (1 - i 



(99) 



which is unity at the associated node and vanishes on 
each face of the cell. Then the shape functions associated 
with the corners of the cell arc of the form 

Niiia) = i [(1 + 6) (1 + 6) (1 + 6) - N,iU] (100) 

where this particular function is associated with the node 
at (1,1,1). The shape functions associated with the other 
corners are generated by the appropriate rotations of this 
function. Again, these shape functions satisfy the criteria 
of locality, partition of unity and continuity. 

Finally, we consider a two-dimensional case that is rel- 
evant for many of the crystal lattices: the square lattice. 



In particular, suppose that the CG region will be treated 
as a two-dimensional projection of the 3D lattice along 
the [001] direction. The square lattice has been used ex- 
tensively in the literature and we include the minimal 
interpolation functions here for reference: 



^i(6) = i(i + 6)(i + 6) 



(101) 



These interpolation functions arc not only useful for 
two-dimensional projections of the lattices treated above 
(simple cubic, fee and bcc with [001] projection), but also 
other more complicated lattices such as the diamond cu- 
bic lattice with [001] projection. 



C. Shape Functions in Reciprocal Space 

In order to make use of the reciprocal space representa- 
tion of the dynamical matrix, it is necessary to have the 
Fourier transform of the shape functions. The Fourier 
transform can be computed numerically, of course, us- 
ing fast Fourier transform (FFT) techniques. In some 
cases it is also possible to compute the Fourier transform 
analytically. In this section, we derive the atomic-index 
Fourier transform of the linear interpolation function in 
one dimension. The result is immediately applicable to 
the 4-node square and the 8-node brick in two and three 
dimensions, respectively. The result is given in Eq. H110() 
below, and readers who are not interested in the deriva- 
tion are free to skip to Section IVl Dl 

Consider the symmetric linear interpolation function 
on a regular mesh in one dimension: 



Xj + l—Xj 



for \x — Xj \ < Xj+i 
otherwise. 



Let a be the lattice constant, and Np, 
We first note the useful identities: 



(102) 
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^ikaNp„/2 sin [ka{Nper + ^)/2} ;^QO-) 



sin [ka{2Nj 



per 



sin (ka/2) 
M)/2] 



sin {ka/2) 



(104) 



that follow from the well known formula to sum geometric 
series [1 + z + + . . . + = (1 - z^+i)/(l - z)] together 
with de Moivre's formula. 

We first transform the atomic index fj, of the shape 
function to the Fourier conjugate variable k: 



(105) 



The shape function Nj^ is expressed as the sum of two 
terms in Eq. (|102|l : the transformation of the first term, 
just equal to unity, is given by Eq. H104|) . but the trans- 
formation of the second is more involved. It is calculated 
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as follows: 
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We have developed an alternative resolution of the zero 
mode problem, which gives an exact formula for the stiff- 
ness matrix with a well defined double inverse (|36|l . Let 
(i'a)/j bc(-'t9fe) /ith component of the ath zero mode of 
D^„; i.e. DtJ.i^i'^a)^' = for a = 1,2,3. Define the 

l^j^fMmode matrices as 
(107) 

(112) 
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Combining the two contributions wc find 
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where iVpe^ is the number of lattice sites per CG cell. In 
higher dimensions, Nper would be replaced by N^^n ^per 
and Np^j.. This result apphes to the regular CG lattice. 
The corresponding formula for a general one-dimensional 
CG lattice is much more complicated, and using a numer- 
ical FFT to calculate it is generally recommended. 



D. The Center-of-Mass Mode 

The stiffness matrix definition involves two matrix in- 
verses. This is somewhat ill-defined because Dfj^i, is sin- 
gular, due to the zero modes. The zero modes are the 
zero energy phonons at the F point in reciprocal space 
that are associated with translation invariance of the cen- 
ter of mass. There are d zero energy phonons in any d- 
dimensional system corresponding to uniform translation 
in each of the d directions. These zero modes make the 
matrix singular. Since we have imposed the criterion that 
the center-of-mass mode should be represented within the 
set of interpolation functions, the singularity is superfi- 
cial. There are two inverses in Eq. H36() . the matrix Kij 
is finite after a suitable regularization. Indeed, the alter- 
nate derivation of the stiffness matrix given in Appendix 
IXI is free from any zero mode problems, so it must be 
possible to devise a suitable regularization scheme. An 
obvious example is 



Kjk = 



1^ I J2 •^J'^ ^^t"" 



ku 
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where /^j^ is the identity matrix. 

The regularization Hlll() is conceptually simple, but in 
practice a small but finite e must be used, and error is 
introduced into the contribution of the long wavelength 
modes. The error can be controlled through the choice of 
an e which is small enough that frequencies of interest are 
not affected appreciably, but large enough that the ma- 
trix is numerically well conditioned. This regularization 
is cheap and adequate for many purposes. 



Using these matrices we construct the projected shape 
matrix 



Then the stiffness matrix 
equation 



Y^V.uNu^ (114) 

fc 

5jk - (115) 

(|36ll is given by the matrix 



(116) 

The nonzero numbers c and c' are arbitrary, but should 
be comparable to the eigenvalues of D^^, to make the 
matrices well-conditioned. This formula works by shift- 
ing the zero eigenvalues in the atomic space and those in 
the nodal space by c and c', respectively, without affect- 
ing the other eigenvalues. The projection matrices undo 
this shift. They are needed within the brackets to stifle 
the crossterms between the zero modes and the nonzero 
modes for incommensurate meshes. For commensurate 
meshes, this formula simplifies to 



commensurate 



X 



= (iViV^) [X -X 
= N{D + c£) ^ 



(118) 



where X • e' is a symmetric matrix since Nif^ (wa)/^ are 
eigenvectors of Xij. 

The zero modes are not integrated out, so a short 
ranged D^i, results in a short ranged Kij. On the other 
hand, a nearest neighbor D^i, does not generally produce 
a nearest neighbor Kij , except where the mesh is atomic 
sized. The stiffness matrix elements typically decrease 
exponentially with separation, so the effective interac- 
tion is short ranged but not nearest neighbor. This is an 
important point, since it is this quality that improves the 
CGMD phonon spectrum. 



VII. (QUASI-)HARMONIC CRYSTALS 

Various properties of harmonic crystals have been com- 
puted within CGMD as a validation of the methodology. 
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A. Phonon Spectra 

The CGMD phonon spectrum offers a good first test 
of the model. Consider a regular, but not necessarily 
commensurate CG mesh. The equations of motion for 
the Hamiltonian (|29() are 

(119) 



where ii" is the nodal acceleration of the j node in the 
a*^ direction. Substitution of a plane-wave normal mode 
M°(t) = yge*''"j ~"^* produces the secular equation 

M(k)c^2^afc = i^°'(k) (120) 

where M(k) and K°''^{k) are the Fourier transform of the 
mass and stiffness matrices, respectively. The form of 
the mass matrix for a monatomic lattice allows further 
simplification: 

m^c^^ ^ [ND'^N'^yJ] (k) (121) 

where we have used Eqs. H32f) and (|37|l . For incommen- 
surate meshes, we have calculated the right-hand side of 
Eq. ((nT|l in real space, then found its Fourier transform 
and solved the secular equation for the phonon frequen- 
cies. One such spectrum is plotted in Fig. 1 of Ref. 17. In 
principle this could be done in any case, but the compu- 
tational cost limits the size of systems that can be treated 
in this manner. We have done calculations with billions 
of atoms per CG cell, but in order to do this it is nec- 
essary to eliminate the real-space representation of the 
dynamical matrix. 

For commensurate meshes with uniform mesh size, we 
may go to reciprocal space and the formulas simplify con- 
siderably. It is even possible to derive analytic formulas 
for the spectra in some cases. Suppose the CG mesh con- 
tains iV^oj^e nodes in the a*'' dimension for a total length 
of La- The CG shape functions in reciprocal space may 
be expressed in terms of a Bravais lattice character x for 
the CG mesh and a CG element structure factor S: 

A^(k,k') = ^e'''-^^-'*''-''^ 

3,1^ 



N. 



(122) 
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where Vlj is any one element from the CG mesh. Note 
that in the atomic limit S is just a delta function in 
the first Brillouin zone. The CG phonon spectrum for a 
monatomic lattice is given by 



i-^|x(k-k')|' |5(k')f 



^'|x(k-k')|' \Sik')\'[D(k')]-' 
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where D{k') is a 3n x 3n matrix where n is the num- 
ber of atoms in the unit cell, and the two inverses are 
matrix inverses. The frequencies are the eigenvalues of 
the resulting matrix. As in Eq. 1)1211) , the first term rep- 
resents the mass matrix in reciprocal space divided by 
rn^; the second term is the middle factor of the stiffness 
matrix. In the atomic limit, the formula reduces to the 
usual expression, D{k)/m. 



B. Analytic Formula for CG Spectrum 

The CGMD spectrum 1)126(1 may be computed in closed 
form for a monatomic solid with a commensurate CG 
mesh in one dimension. We presented the analytic ex- 
pression for the spectrum with nearest neighbor interac- 
tions and linear interpolation in Ref. IT^ Eq. (12): 
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where the sums over p run from to Np 
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1, Nr, 



^atomZ-^nodc IS the number of atoms per cell and K is the 
nearest- neighbor spring constant. This formula shows 
the contribution of many modes of the underlying crys- 
tal to each CGMD mode, resulting from the choice of 
interpolation functions which have many normal mode 
components. Near the center of the CG Brillouin zone, 
a single mode (p = 0) dominates the sums (|127|) . This 
dominance reflects the fact that long wavelength modes 
are well represented on the CG mesh. Near the bound- 
ary of the CG zone [k « N^ode''^ /{Na)], many modes 
contribute. The many modes are needed because peri- 
odicity forces the slope of the spectrum to zero, and the 
modes act in concert to keep the CGMD spectrum close 
to the true spectrum which is not smooth at the bound- 
ary. For comparison, the formula for the lumped mass 
FEM spectrum is 
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the formula for the FEM spectrum with the distributed 
(consistent) mass matrix is 



OJ 



dist 



(fc) 




sin I ^kNperC 



and the formula for the exact MD spectrum is 



(129) 



,MD 



(k) = 2W — sin ifc. 



(130) 



The coarse-grained mass and stiffness matrices conspire 
to produce a Fade approximant of the true spectrum, 
and thereby achieve the O(fc^) improved relative error, 
compared to the 0{k'^) relative error of the two FEM 
spectra. 
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The remainder of this subsection is devoted to the 
derivation of the analytic formula for the CGMD spec- 
trum H127|l . We make use of the formula for the Fourier 
transform of the symmetric linear interpolation function 
on a regular mesh in one dimension, given above in Eq. 

(irrnii : 



N,{k') 



sin^(fc^QiVper/2) ,fc'^. 



Nper sin^(fc'a/2) 



where in this formula only the atomic index has been 
transformed. The Fourier transform of the index j is 
straightforward, and the spectrum could then be derived 
using Eq. (|126|l . We take a different approach. The spec- 
trum is given by 
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M{k) 

1 / {NNT)^ 




[ND-^NT\ 



(131) 
(132) 



where we have made use of the formulas for the stiff- 
ness and mass matrices, Eqs. (|57|) and (15^ . respec- 
tively. The subscript k denotes the Fourier transform 
of wavenumber k, and we note that the Fourier trans- 
forms of both the nodal indices and the atomic indices 
take on values of the form 27rn/i, but for the atomic 
indices —^L/a < n < ^L/a whereas the Fourier trans- 
form of the nodal index lives in a reduced Brillouin zone, 
— ^^nodc < n < ^iVnodc: whcrc iVnodc IS the number of 
nodes (in one dimension). 

Evidently, we need to calculate quantities of the form 
(A^XiV^) ^, where the matrix X is either the identity ma- 
trix or the inverse of the dynamical matrix. Such quan- 
tities are calculated in the following way: 

{NXN^)^ = 7V-;,^,^e-'=(---^)(7VXArT),^/133) 
hi 
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where Eq. (|134|1 follows from Eq. (|ll()|l , and we have used 
Xi — Xj — {Aj)Nper-a where Aj = To get from H134|) 
to H135() , we have used the fact that the sum over Aj gives 
a delta function in the reduced Brillouin zone; i.e., a sum 
of delta functions periodically repeated through the full 
Brillouin zone. 



The CGMD spectrum is then calculated using Eq. 
(|132|l together with Eq. I|135|l with X equal to the iden- 
tity in the numerator and D^^ in the denominator. The 
result is 
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where D{k) is the dynamical matrix in fc-space and the 
sums over p run from to Nper- For a nearest-neighbor 
harmonic model D{k) = AK sin^ {^ka) . Substitution of 
this into Eq. H13t)|) results in the analytic formula for 
the CGMD spectrum that appears above (|127|l . This 
calculation may be generalized to 3D, where D{k) is a 
3x3 matrix and the shape functions are products of linear 
interpolations functions in each of the three dimensions. 

The first test is the phonon spectrum for atoms with 
harmonic interactions coarse grained to a regular, but not 
necessarily commensurate mesh. The normal modes are 
plane waves both on the underlying ring of atoms and on 
the CG mesh. The wave vector k is a good quantum num- 
ber for both. The nonzero terms of the dynamical matrix 
are of the form: D^^ = 2K, Dp^p±i = —K. Figure 1 of 
Ref.T? shows the resulting phonon spectra in four cases: 
exact, CGMD, distributed mass FEM and lumped mass 
FEMiffi The latter two use the long wavelength elastic 
constants. The spectra are for a periodic chain of 1024 
atoms with lattice constant a coarse grained to 30 nodes. 

Figure 1 of Ref. [l3 shows that CGMD gives a bet- 
ter approximation to the true phonon spectrum than the 
two kinds of FEM do. All three do a good job at the 
longest wavelengths, as expected, but CGMD offers a 
higher order of accuracy. The relative error for CGMD 
is 0{k'^) while that of the two versions of FEM is only 
0{k^). At shorter wavelengths, there are significant devi- 
ations from the exact spectrum. The worst relative error 
of CGMD is about 6%, better by more than a factor 
of three than that for FEM. This improvement is made 
possible by the longer-ranged interactions of CGMD as 
compared to FEM. The continuity condition satisfied by 
linear interpolation is enough to ensure that the hydro- 
dynamic modes {k 0) are well modeled, but the lack 
of continuity of the derivatives shows up as error in the 
spectrum of the modes away from the zone center. This 
error vanishes for the smooth, nonlocal basis consisting of 
the longest wavelength normal modes. It turns out that 
the CGMD error at the CG zone boundary is relatively 
small (less than 1%) for technical reasons. Also note that 
even though the number of atoms varies from cell to cell 
in the incommensurate mesh, the CGMD spectrum is 
free of anomalies. Other computations have shown that 
CGMD with linear interpolation is well behaved on ir- 
regular meshes, as well. 

We are now in a position to investigate the effect of the 
CGMD relaxation terms eliminated in R-CGMD. They 
make the CGMD stiffness matrix non-local and therefore 
add to the cost of CGMD. What is the benefit of this ad- 
ditional computational complexity? Using the same pro- 
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FIG. 5: A comparison of the acoustic wave spectra for CGMD, 
R-CGMD and MD with frequency vs. wavenumber plotted. 
The units are ko = ^/{Npera.) and ujo = fcofl-y/ K/m. The 
MD spectrum, corresponding to a simple ID ball and spring 
model, is the ideal case. The CGMD and R-CGMD spec- 
tra are computed on a regular mesh with Nper = 32 atoms 
per cell. Both full CGMD and its rigid approximation, R- 
CGMD, are in good agreement with the MD spectrum for 
fc ~ 0; i.e., at long wavelengths. At short wavelengths near 
the zone boundary, the CGMD spectrum is more accurate 
than the R-CGMD spectrum, a property that we attribute to 
the relaxation effects accounted for by CGMD but eliminated 
in the rigid approximation. 



cedure outlined above, we have computed the R-CGMD 
spectrum: 
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where the two lines are to be compared with the CGMD 
results (fna and ifT^ . respectively. The CGMD and R- 
CGMD spectra are plotted together with the MD spec- 
trum in Fig.[Sl It is clear that both approaches work well 
for long wavelengths (fc ^ 0). In fact, it is reasonable that 
the relaxation should be unimportant in this regime since 
the displacement field is varying slowly on the scale of the 
mesh, so the lowest energy configurations of the MD best 
fits to the interpolated displacement field should be close 
to having the atoms at their linearly interpolated posi- 
tions. At short wavelengths (near the CG zone bound- 
ary), the story is different, however. The displacement 
field is varying at the scale of the mesh, and the atoms can 
reduce the energy through relaxation. This is evident in 
the improved value of the zone boundary frequency for 
CGMD (0.67% error) vs. R-CGMD (10.2% error). As 
shown in Fig. (HJ these values are typical, and close to 
the asymptotic value for coarse meshes. This difference 
in the performance of CGMD and R-CGMD is evident 
in many properties sensitive to the coarse-grained lattice 
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FIG. 6: A comparison of the error in the acoustic wave fre- 
quency at the zone boundary [k = t: / {Npera)] for CGMD 
and R-CGMD as a function of the level of coarse graining, as 
expressed by Np^r- At Nper ~ 1, there is no error in either 
frequency. At all other values of Np^r, the error in the CGMD 
frequency is less than that of R-CGMD, asymptotically going 
to 0.66% and 10.3%, respectively. 



dynamics at the zone boundary, such as the scattering 
properties we consider below. It is interesting to note 
that while the magnitude of the error is appreciably dif- 
ferent, it is small in both cases. This suggests that there 
may be ways to formulate an approximation to CGMD 
that is intermediate between CGMD and R-CGMD, both 
in terms of the detail of the physics that is described and 
the computational cost. This is a topic we will address 
in the future. 



C. Numerical Calculation of Argon CG Spectrum 

It was shown in Ref. ^3 that the CGMD phonon spec- 
trum is closer to the true spectrum than that of FEM for 
a one dimensional chain of atoms with nearest-neighbor 
interactions. We now compare phonon spectra for a three 
dimensional real material: solid argon. Some of these re- 
sults were reviewed briefly in Ref. Ij. We also treat tanta- 
lum below. We find that CGMD again offers an improved 
spectrum. 

Solid argon crystallizes in the fee structure where it is 
well-described by the Lennard-Jones potential 
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where e = 1.63 x 10"^^ J and cr = 3.44 A.*^^ The 
elastic constants for this potential are given by Cn = 
105.3e/cr3 = 4.21 GPa and C12 = C44 = 60.18e/cr''' = 
2.41 GPa. 

We use the four-fold symmetrized 8-node brick inter- 
polation functions defined on the fee Bravais lattice for 
both CGMD and FEM. The fee lattice is generated by 
the unit cell vectors ai = (0,§,|), a2 = (|,0,§), and 
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aa = (|,§,0). As discussed above, the conventional 
rhombohedral interpolation functions of the 8-node brick 
consist of products of one- dimensional linear interpola- 
tion functions: 



7V,(x) = n^[2ba 



(140) 
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as in Eq. H94() These interpolation functions break the cu- 
bic symmetry. Even though the effect of the symmetry 
breaking is small, it complicates the analysis of the spec- 
trum, causing deviations from the true spectrum that 
have nothing to do with the intrinsic accuracy of the 
methods. 

In order to eliminate the small symmetry-breaking ef- 
fects completely, we introduce a symmetrized version of 
the interpolation functions that is useful for spectrum 
calculations. We construct a multiplet of interpolated 
fields, where each component of the multiplet is an image 
of the mesh associated with Eq. (|140|l under the action 
of the point group. For the fee lattice with rhombohe- 
dral cells, there are four inequivalent transformations of 
the mesh, so our field becomes a four component vector. 
The frequency values are then averaged over these four 
components. This is equivalent to the standard group 
theoretic operation of averaging over the orbit in order 
to restore symmetry. The four inequivalent group oper- 
ations are the C2 elements 



±10 
3=1 ±1 
±1 



(141) 



with —1 appearing an even number of times such that 
5 is a proper rotation with det(g) = 1. Another way 
to view this symmetrization procedure is that a sym- 
metrized stiffness matrix is used 
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where the sum is over the four diagonal matrices g (|141|l . 
If the meshing did not break the point group symmetry, 
there would be no need for symmetrization and, indeed, 
the sum in Eq. (|142|l would reduce to a single term. The 
symmetrization procedure is simply a way to restore the 
symmetry broken by the mesh in order to facilitate anal- 
ysis and comparison of the spectra. 

We now undertake the actual calculation of the stiff- 
ness and mass matrices. The FEM matrices are com- 
puted as follows. The mass matrix takes one of two forms. 
The mass matrix computed strictly from the interpola- 
tion functions is known as the distributed or consistent 



mass matrix. It is given by 



M, 



dist 



cPxpNi{x)Nj{x) 



(143) 



d^N,{^a)N,{^,) (144) 



miV, 



per 



2^3 I I'^-^i 

6 



(145) 



where / (0 < / < 3) is the number of edges in the shortest 
path along a single hexahedral cell connecting nodes i 
and j. The number Np^r is a generalization of the one- 
dimensional case, where now it is the number of lattice 
sites along one dimension of the cell, and for cells with 
different dimensions in the three directions N^^^ should 
be replaced by N^^^N^^^Np^j.. The Fourier transform 
of M!^j'^^ may then be calculated as a sum over the 27 
neighboring nodes (indexed by ni, 712,713 running from 
-1 to 1) 
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where a(, is the 6*'* real space basis vector. 

For many applications, it is sufficient (and in some 
cases even more accurate) to use a diagonal approxima- 
tion to the mass matrix known as the lumped mass ma- 
trix. It is given by 
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so each element of the diagonal is just equal to the mass 
contained in the Voronoi cell about the corresponding 
node. The Fourier transform of M'"™^ is 



M'""^P(k) 



(148) 



Note that M'^'"'(k = 0) = M'"'"P(k = 0), and they are 
equal to the Voronoi mass as they should be. 
The FEM stiffness matrix is given by 
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where the prefactor a is the lattice constant of the rhom- 
bohedral unit cell and the elastic tensor Caked has been 
contracted with the reciprocal space metric as appropri- 
ate for the non-orthogonal coordinates. This equation 
has too many components to present the complete ex- 
pression here [cf. Ref . Is^ . Nevertheless, the calculation 
is elementary algebra, and the results were used to cal- 
culate the Fourier transform. The result is a 3 x 3 matrix 
for each value of k: Kj^^^'^ {k). 
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The formulas for the CGMD mass and stiffness matri- 
ces in real space for rhombohedral elements are computed 
similarly. The mass matrix is given by 
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found using Eq. H110() to be 



b=ipt^i^P- sin* (afe -1^/2) 
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where Li, is the length of the CG cell in the 6*'' direction. 
Upon substitution back into Eq. H153(l . we find that 



where as in Eq. p45l) . ^ (0 < Z < 3) is the number of 
edges in the shortest path along a single rhombohedral 
cell connecting nodes i and j. The correspondence of 
the leading terms to the terms in the FEM distributed 
mass matrix (|145|) is evident, so that CGMD reproduces 
the FEM distributed mass matrix in the large- iVper limit. 
This expression assumes that the mesh consists of trigo- 
nal cells in which the linear dimensions are equal in all 
three dimensions, but it could be generalized immedi- 
ately to unequal dimensions. The Fourier transform is 
given by 
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where as in Eq. H14t)|l . the 27 neighboring nodes are in- 
dexed by n\,n2,n-i running from -1 to 1. Note that the 
CGMD mass matrix also satisfies the mass sum rule: 
MCGMD(j^^0) = m7V3^,. 

The CGMD stiffness matrix is calculated according 
to Eq. H37|l using the reciprocal space representation 
of the dynamical matrix. In particular, we calculate 
the spectrum using Eq. H132() . suitably generalized to a 
monatomic lattice in three dimensions: 

u:U^) = ^M(k) {ND-^N^)-\V) (153) 

where the actual frequency on each phonon branch is 
given by the square root of one of the three eigenvalues 
of the 3x3 matrix ^^^^(k). Here we have made use of the 
expression for the mass matrix of a monatomic lattice 
()32|l . The denominator is part of the stiffness matrix. 
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where D'^^i^') is the matrix inverse of the 3x3 matrix 
-Dafc(k'). The outer inverse is a 3 x 3 matrix inverse, as 
well. The Fourier transform of the shape functions is 
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where kp = k + 2'!ipi)/Li). The mass matrix Af (k) is 
given by Eq. (|152|l . It is in the numerator, what might 
seem to be the wrong place, because of the form of the 
monatomic secular equation (|121|l . It is clear from Eq. 
(|157|l that CGMD reproduces the MD spectrum in the 
long wavelength limit 
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which follows from expanding sin(j;) = x -\- . . . for small 
arguments. In the short wavelength limit the many terms 
in the sum over pf, contribute to Eqs. H15t)|) and (|157|l . 
ensuring periodicity in the CG reciprocal space. 

The spectrum is then computed in each case from the 
resulting secular equation at each value of k. For the true 
spectrum, the secular equation is 

det [w2(k) - Dab{\i)/m\ - 0, (159) 

for the CGMD spectrum the secular equation is 



det [w2(k)^„,_c^2^(k)] =0 



(160) 



with w^j, given by Eq. H157|) . and for the FEM spectra 
the secular equation is 

det J^, _ X,b(k)/M(k)] = (161) 

with M (k) given by Eqs. ifT^ and for distributed 

and lumped mass, respectively. The determinant of the 
3x3 matrix gives a cubic secular equation with three 
eigenvalues. The eigenvalues are real for all of the cases 
considered here. The Lennard- Jones potential for argon 
was cut off in real space after the twelfth nearest neighbor 
shell for both the MD and CGMD spectra. For many ap- 
plications it would not be necessary to extend the range 
this far; however, in this case we wanted to test a po- 
tential extending well beyond nearest neighbors. The 
semi-analytic formulas that we have presented here, such 
as Eq. H157I) . have been used in order to calculate the 
spectra of extremely large systems with minimal com- 
putational expense, and all of the plots presented be- 
low were calculated using these equations. It should be 
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FIG. 7: The phonon spectra for solid argon in the atomic limit 
are shown as plots of wave frequencies vs. wave number along 
various high symmetry lines through the Brillouin zonei^ The 
high symmetry points are labeled according to the standard 
convention for an fee latticeSS*^; for example, F is at the cen- 
ter of the zone, k = and X = (l,0,0)7r/a where a is the 
fee lattice constant. The 9 curves represent the 3 branches 
of the spectra for CGMD and two versions of FEM, one us- 
ing the conventional distributed mass matrix and one using 
a diagonal lumped mass matrix. In the atomic limit CGMD 
reproduces the Lennard- Jones MD spectra exactly, whereas 
the FEM spectra show significant error, especially with the 
distributed mass matrix. 

emphasized, however, that the spectra could have been 
calculated using the basic real-space matrix formulas or 
using numerical Fourier transforms. This has been done 
for the smaller systems as a validation. 

The spectra have been computed for three levels of 
coarse- graining: the atomic limit (1x1x1 or no coarse- 
graining), a shght coarse-graining (2x2x2) and a case 
approaching the continuum limit (32x32x32). The num- 
bers Ux X Uy X Uz indicate the number of atoms within 
an unsymmetrized CG cell in each direction, i.e., Nper = 
1, 2, 32, respectively. These values correspond to cells 
containing 1, 8 and 32678 atoms, respectively. The re- 
sults are shown in Figs.[7||Hland|Hl 

Elastic wave spectra are of interest because they pro- 
vide a means of quantifying the small-amplitude dynam- 
ics of the system. They represent the energetics of every 
normal mode of vibration of a system of atoms. In coarse- 
grained dynamics, the wave spectra provide an excellent 
way to quantify the fidelity of the coarse-grained model. 
Since the normal modes of a crystal are plane waves, they 
are uniquely identified by a wave number k and a branch 
index, for example indicating a transverse optical mode 
or a longitudinal acoustic mode. The normal modes cor- 
respond to a lattice of points in reciprocal space (k-space) 
inside a bounded region known as the Brillouin zone. It 
is not possible to plot ci;(k) for k throughout the three 
dimensional Brillouin zone, so typically a'(k) is plotted 
along lines, in particular high symmetry lines, through 
the 3D Brillouin zone.^M2, This convention has been used 



FIG. 8: The phonon spectra for solid argon on a mesh with 
slight coarse graining are shown as plots of wave frequencies 
vs. wave number along various high symmetry lines through 
the Brillouin zone (cf. the caption to Fig. |7J. The 12 curves 
represent the 3 branches of the spectra for MD, CGMD and 
two versions of FEM, one using the conventional distributed 
mass matrix and one using a diagonal lumped mass matrix. 
Each cell of the CG mesh contains 8 atoms. For this first level 
of coarse graining the CGMD spectra is in better agreement 
with the MD spectra than either of the FEM spectra. Of 
the two FEM spectra, the lumped mass spectra is somewhat 
better. 



in Figs. El Hand El 

Consider the spectra in the atomic limit shown in Fig. 
[3 The CGMD spectrum agrees precisely with the true 
spectrum. Of the two FEM spectra, the lumped mass 
spectrum is closer to the true spectrum. This is sen- 
sible because the mass is localized to the nodes in the 
atomic limit, since each node represents one atom. Over- 
all, the lumped mass frequencies are lower than the true 
frequencies, whereas the distributed mass frequencies are 
higher. This ordering remains true regardless of the level 
of coarse-graining. 

In the continuum limit shown in Fig. El the CGMD 
spectrum no longer agrees exactly with the true spec- 
trum, but it is still a better approximation than either 
FEM spectrum. It is clear that in the continuum limit, 
the distributed mass produces the better spectrum of the 
two finite element cases. This again makes sense, because 
the mass is becoming more evenly spread throughout the 
CG cell. Still, the CGMD spectrum is significantly better 
than the distributed mass FEM spectrum. 

The intermediate case is shown in Fig. |S1 Already the 
distributed mass FEM is in better agreement with the 
MD spectra than the lumped mass FEM is. It is remark- 
able that cell containing as few as 8 atoms are beginning 
to exhibit continuum behavior. This one-two-many qual- 
itative dependence is typical of many large- A'^ expansions, 
where the large- TV limit quickly becomes a good approx- 
imation to the real system behavior, and even for TV as 
low as two or three it is a good approximation. The 
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FIG. 9: The phonon spectra for solid argon on a mesh ap- 
proaching the continuum limit (32768 atoms per cell) are 
shown as plots of wave frequencies vs. wave number along 
various high symmetry lines through the Brillouin zone"' (cf. 
the caption to Fig.[7||. The 12 curves represent the 3 branches 
of the spectra for MD, CGMD and two versions of FEM, one 
using the conventional distributed mass matrix and one using 
a diagonal lumped mass matrix. With this significant level of 
coarse graining the CGMD spectra is again in better agree- 
ment with the MD spectra than either of the FEM spectra. 
Of the two FEM spectra, now the distributed mass spectra is 
better. 



CGMD spectra are again in better agreement with the 
MD spectra than either of the FEM spectra are. 

It should be emphasized that the shape of the FEM 
spectra is the same in the three plots. Continuum elas- 
ticity is scale invariant, and the changes in in FEM spec- 
tra are a simple rescaling of frequency and wave num- 
ber. This scaling is clearly evident in Eq. (|146|l for the 
FEM distributed mass. The scale, NpeT^ enters through 
the prefactor TV^g^ scaling the frequency and the factor 
of A^per in the argument of the cosine scaling the wave 
number. The same scaling of the wave number is present 
in the stiffness matrix, but its prefactor goes like Nper 
rather than A^pg^- a result, the frequency (~ ^JK/M) 
and the wave number scale like l/Nper- In the linear 
portion of the spectra near fc = 0, the two effects cancel, 
but the spectra are modified significantly near the zone 
boundary. Thus, when we discuss the comparison of the 
MD spectra with the FEM spectra and find better agree- 
ment with the lumped mass FEM for small cells and bet- 
ter agreement with distributed mass FEM for large cells, 
it is not that the FEM spectra are changing. The MD 
spectra are changing from discrete atomic behavior to 
continuum behavior. It is the natural scale dependence 
of true lattice dynamics. CGMD has this scale depen- 
dence arise naturally, as well, and so it is able to a large 
extent to track the changes in the MD spectra with cell 
size. 

If we look at the differences between the spectra in 
more detail, we can start to understand how the under- 



lying physics gives rise to these differences. On all three 
plots, the transverse phonons are degenerate (have the 
same frequency) along the lines connecting the zone cen- 
ter and the middle of the zone faces, F — X in the [100] 
direction and F — i in the [111] directioniS This degen- 
eracy follows from the C4 and C3 symmetry along those 
lines, respectively. The lines along the zone boundary 
and the T — K line have reduced symmetry (at most 
C2), and the transverse phonons are not degenerate in 
those cases. There are some isolated points of increased 
symmetry. The W point is an example. One transverse 
phonon and the longitudinal phonon are degenerate at 
W = (2, 1,0)^, as can be seen in the MD curve in Fig. 
This happens because W is the mid-point on a line 
between X on two adjacent Brillion zones, and the lon- 
gitudinal mode at one X becomes the transverse mode 
at the other X, and vice- verse. The two branches must 
cross at the mid-point, and hence they are degenerate at 
W . This is no longer true for the MD phonons in the 
coarse grained systems, since Wcg no longer has this 
mid-point property. As a result, the frequency of the 
longitudinal MD mode increases going from LtoW {W 
is farther from F), whereas it decreases for CGMD and 
the two versions of FEM since the longitudinal frequency 
is dropping to meet the transverse frequency. This is the 
most pronounced discrepancy in the spectra. It is par- 
ticularly bad for the lumped mass FEM as the cell be- 
comes large. CGMD is in better agreement with the two 
transverse MD modes at W , while the distributed mass 
FEM, which tends to be high overall, is in better agree- 
ment with the high longitudinal MD frequency at W . In 
general, the CGMD spectra are in substantially better 
agreement with the MD spectra than either of the FEM 
spectra are. 

It is not obvious, but even in the long wavelength limit 
(near the F point), the CGMD spectrum is better than 
either of the two FEM spectra. The relative error for the 
CGMD spectrum is of order 0(fc''), whereas it is O(fc^) 
for the two FEM cases. This improved error was demon- 
strated above with the formulas for the ID frequencies, 
and it continues to hold in 3D. An order 0{k^) relative 
error is the naive expectation, since the phonon disper- 
sion relation is linear, with the leading corrections of or- 
der 0{k^) due to symmetry. The higher order error for 
CGMD is due to a subtle cancellation between the mass 
and stiffness matrices. This cancellation can be seen from 
the formula for the general ID CGMD spectra ()136|l and 
its 3D generalization 

^2 (k) ^ lyr Ep?risin-^(la,.kpJ 

where kpj_ = k + -^r\^b and b;, is the reciprocal lat- 
tice basis element. Note that this formula is equivalent 
to Eq. H157() . which is the equation we actually used to 
compute the CGMD spectra. The two equations differ 
because they make use of different expressions for the 
mass matrix. 
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It is interesting again to compare the spectra for 
CGMD and R-CGMD, now for the 3D phonons. In the 
atomic hmit {N^^^ = 1); the two agree with each other 
and with the exact MD spectrum. For coarsened lat- 
tices (N^;^^ > 1), near the T point (k = 0), CGMD and 
R-CGMD are in good agreement (R-CGMD not pfotted 
here), as was observed in the ID case and shown in Fig. 
|S1 We may compare the formulas for the spectra, Eq. 
(ITH^ and its R-CGMD counterpart, 




The two equations for the spectra may be expanded 
about the T point, and both exhibit the improved rel- 
ative error, 0{k'^). Near the zone boundary, R-CGMD 
behaves more like the distributed mass FEM case. It is 
here that the relaxation physics built into CGMD has its 
most pronounced effect, especially at the high symmetry 
point L on the acoustic branch, where the full CGMD 
error is very small. 



D. The Finite Temperature Tantalum CG 
Spectrum 
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FIG. 10: The room temperature phonon spectra for solid tan- 
talum on a mesh with no coarse graining are shown as plots 
of wave frequencies vs. wave number along various high sym- 
metry lines through the Brillouin zone (for comparison to the 
solid argon case shown in Fig. . As in the case of argon, 
the CGMD phonon spectrum agrees exactly with the MD 
spectrum when the mesh is refined to the atomic limit. It is 
common practice to leave the gap between the second N and P 
points since that part of the spectrum is already represented 
to the left. 



We have calculated the elastic wave spectrum for a va- 
riety of materials. As a second example we present results 
for the phonon spectra of tantalum at room temperature. 
Tantalum was chosen to demonstrate CGMD in a more 
open crystal structure (bcc) and in a system using many- 
body interatomic potentials. We use the Finnis-Sinclair 
many-body potential for tantalum^ with the improved 
Ackland-Thetford core repulsiouiS^ The elastic constants 
for this potential are Cn = 266.0 GPa, C12 = 161.2 GPa, 
and C44 = 82.4 GPa. We calculate the finite temperature 
dynamical matrix in a conventional MD simulation con- 
sisting of 2000 atoms in a lattice of 10x10x10 bcc unit 
cells with periodic boundary conditions. The system is 
equilibrated to T = 300 ± O.IK and P = ± 10"^ Qp^ 
through scaling of the box size and velocities every 100 
time steps until the target temperature and pressure 
were attained and then an additional 5000 steps without 
rescaling to ensure equilibration. The equilibrium lattice 
constant at this temperature was found to be 3.3129 A, 
expanded by 0.2% from the T = OK value of 3.3058 A. 

Subsequently, the dynamical matrix was calculated ev- 
ery 1000 time steps averaged over every atom in the sim- 
ulation. With the Finnis-Sinclair potential it is possible 
to use an analytic expression for the dynamical matrix, 
since it is possible to take two derivatives of the energy 
analytically with respect to atomic displacements. The 
expression is given in Ref . l64l In principle we are com- 
puting an ensemble average, which we have implemented 
by averaging over the equivalent lattice sites of the sys- 
tem and over multiple time steps (relying on ergodicity 
for the equivalence of ensemble and temporal averages 
in equilibrium). In all, we have averaged over a total of 



10 snapshots of the system and imposed the cubic {Oh) 
point group symmetry by averaging over the point group 
operations. The range of the dynamical matrix includes 
out to the sixth nearest neighbor shell in tantalum at 
T=300K (the range of the pairwise functions entering the 
potential includes the first and second nearest neighbor 
shells) . 

The results for tantalum are very similar to those for 
solid argon in terms of quality. The elastic wave spectra 
are plotted in Figs. 1101 andllll The first figure shows the 
spectrum of elastic waves on a fully refined mesh. The 
CGMD and MD spectra agree exactly and are overlap- 
ping. The second figure shows the spectrum on a mesh 
consisting of 8 atoms per rhombohedral cell, as described 
in the argon case. Here the spectra do not agree exactly, 
but the results are comparable in quality to those from 
the argon simulations. In comparing to the correspond- 
ing argon plots (Figs. [71 and (S)) it should be noted that 
the high symmetry k points are somewhat different due 
to the differences in the bcc and fee crystal structures. 
For example, the F-H line is in the [100] direction and the 
F-P line is in the [111] direction for bcc Tantalum. Again 
CGMD performs better than either of the two FEM mod- 
els (data not shown) in reproducing the MD spectra. The 
many-body potential and the more open crystal structure 
do not have a significant impact on the quality of the re- 
sults. 
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FIG. 11: The room temperature phonon spectra for solid tan- 
talum on a mesh with 8 atoms per mesh cell are shown as plots 
of wave frequencies vs. wave number along various high sym- 
metry lines through the Brillouin zone (cf. Fig.lHJ. As in the 
case of argon, the CGMD phonon spectrum agrees very well 
with the MD spectrum in the limit of long wavelengths (near 
the r point), and it agrees reasonably well throughout the 
coarse-grained Brillouin zone. 



E. Dynamics and Scattering 

Most of the applications that we envision for CGMD 
are dynamical and have varying mesh size. For example, 
in studies of crack propagation it may be advantageous 
to introduce a coarse-grained model of the far-field re- 
gions away from the cracki^^ For these applications it is 
important that elastic waves generated at the crack tip 
do not reflect from the coarse-grained region and per- 
turb the crack propagation,^&SLSigi2a CGMD offers the 
chance to allow the longest wavelength modes to propa- 
gate much farther into the periphery within incurring a 
commensurately greater computational expense. Other 
coarse-grained models of the periphery, such as hybrid 
FEM/MD schemes, may offer the same advantage. In 
energy-conserving CGMD, the short wavelength modes 
are reflected from the CG region, but this process is 
sufficiently dispersive that shock waves are smoothed 
out and the potential wave-reflection pathologies are 
averted. The unphysical wave reflection may also pro- 
duce a nonzero Kapitza resistance at the interface, which 
can lead to an unphysical temperature gradient across 
the interface. Of course, a stationary system started in 
thermal equilibrium remains at a constant, uniform tem- 
perature given a reasonable measure of temperature in 
the CG region, but a system driven out of equilibrium 
may exhibit unphysical gradients on time scales shorter 
than the relaxation time. 

Given the potential problems associated with wave re- 
flection, we have developed a methodology to quantify 
the problem. The natural measure of the ability of waves 
to propagate from an atomistic region into a CG region 
is the S'-matrix of scattering theory, or in one dimen- 



sion, the transmission and reflection coefficients, T and 
7?., respectively. The basic approach to scattering prob- 
lems is to look for solutions of the equations of motion 
of the form of an incoming plane wave and an outgoing 
spherical wave. 
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where this asymptotic form of the displacement field 
holds well outside the scattering region. The S'-matrix 
and the scattering cross section may be determined from 
/k(k). For CGMD, the scattering region is the region 
where the stiffness matrix differs from the MD dynami- 
cal matrix. A tremendous amount of theoretical analysis 
has been developed for scattering problemsi^i In lattice 
dynamics, scattering is complicated by the anisotropy of 
the lattice. The asymptotic form given above (|164|l is 
only applicable to isotropic scattering, but the formalism 
is readily generalized. To the best of our knowledge, 3D 
scattering cross sections have not been computed for any 
of the proposed solutions to the wave reflection problem. 

We restrict our discussion to the one-dimensional case, 
for which the analysis is more straight-forward. We 
have calculated these scattering properties for CGMD 
and FEM, for comparison. The reflection coefficients are 
computed in the same way for both. The asymptotic re- 
gion is described by harmonic MD on a regular lattice, 
and the normal modes are the well-known plane wave 
solutions. As in Eq. (|164|l . the asymptotic form of the 
displacements is known for each frequency: 
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) for j < 1 
for j > N 



(165) 

where the reflection and transmission coefficients are 
given by R = |rp and T = respectively. We have 
assumed that the scattering region is contained between 
xi and xn, in the sense that these points bound the 
CG region of the mesh and are separated by more than 
the range of the MD potential from any coarse-grained 
cell. This requirement guarantees that the form of Uj(t) 
in the relation (|165|) is a strict equality and not just an 
asymptotic relation like H164() . The wave number k is 
determined by the frequency, mw^ = D{k), where D{k) 
is the MD dynamical matrix. The leading coefficient A 
just determines the amplitude and is irrelevant for our 
purposes, so we set ^ = 1. In principle, the displace- 
ment field could have components with many different 
frequencies, but since the problem is linear, we may re- 
strict to a single frequency without loss of generality. 
Note that while we are considering the harmonic prob- 
lem with a short-range MD potential, this analysis could 
be generalized to non-linear wave propagation and non- 
linear or long-range scattering using the LSZ scattering 
formalismf£^ The incoming and outgoing waves forming 
the asymptotic boundary conditions (|165|l would need to 
be suitably dressed. Then the scattering cross section 
could be expressed in terms of the one-point irreducible 
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Feynman graphs. This extension could be interesting, 
but the Hnear case will suffice for our purposes. 
The equations of motion are given by 

M,jUj{t) = -Kauiit) (166) 

and we look for solutions with angular frequency lo, 

such that the asymptotic boundary conditions (|165|l are 
satisfied: 



r e for j < 1 

(t-l) ^ for j > N 



(168) 



where again these boundary conditions are a strict equal- 
ity. The equations of motion for Vj are 



(169) 



In principle, there are many ways to solve the equations 
of motion 1)169(1 with the boundary conditions 1168(1 : in 
practice, we found this problem to be rather subtle. A 
similar scattering problem must have been solved before, 
but we have not been able to find a solution in the liter- 
ature. The approach we take here is to note that we can 
relate the solution in the outer regions to the solution at 
the boundary points 
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where n > 0. Here a is the lattice constant. Using 
this trick, the problem is reduced to the calculation of 
wi,. . . jfAT using the following N equations: 

{K,,-uj^Nhj)vj = -(A'y-c^2^/,j)e''="^ (172) 
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where in practice the sums just run out to the range of 
the potential. The implicit sums over j run from 1 to 
A''. Then the scattering coefficients are determined by 



7^ = 



and T = \tr with 
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which follow from Eq. I(168() . 

In Fig. El we plot the reflection coefficient TZ{k) for 
scattering from a CG region of 72 nodes representing 
652 atoms in the middle of an infinite harmonic chain 
of atoms. The reflection coefficients for CGMD, lumped 
mass FEM and distributed mass FEM are plotted. The 
lattice is shown in Fig. E| The cell size increases 
smoothly in the CG region, as it should, to a maximum 
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FIG. 12: A comparison of the reflection of elastic waves from 
a CG region in three cases: CGMD and two varieties of FEM. 
Note that the reflection coefficient is plotted on a log scale. 
A similar graph plotted on a linear scale is shown in Ref. 
The dashed line marks the natural cutoff [feo = / {Nmaxa)], 
where Nmax is the number of atoms in the largest cells. The 
bumps in the curves are scattering resonances. Note that at 
long wavelengths CGMD offers significantly suppressed scat- 
tering. 



of Nmax = 20 atoms per cell. In all three cases shown TZ 
becomes exponentially small in the long wavelength limit, 
and it goes to unity as the wavelength becomes smaller 
than the mesh spacing — a coarse mesh cannot support 
short wavelength modes. The threshold occurs approxi- 
mately at = tt/ {Nmaxo), the natural cutoff correspond- 
ing to a wavelength of A = 2NmaxO,- The threshold for 
CGMD takes place almost exactly at this point because 
the CGMD phonon frequencies are more accurate than 
those of the two versions of FEM. According to three di- 
mensional scattering theory in the limit of wavelengths 
much larger than the size of the scattering region, the 
scattering cross section should vary like a ^ k"^. This 
favorable transmission of long wavelengths is the well 
known Rayleigh scattering that gives us a blue sky.^^ In 
these one dimensional scattering calculations, the trend 
is for TZ ^ where the exponent P is roughly /3 « 4 ± 1 
for FEM and /3 « 8 ± 2 for CGMD. We hypothesize that 
the difference is due to the improved accuracy of CGMD 
at long wavelengths. Using the Born approximation, the 
scattering strength should go roughly like r ~ fc''', where 
7 is the order of the error, or 7 = 2 in FEM and 7 = 4 
in CGMD. Then the reflection coefficient would go like 
n = |r|2 - fc27^ giving /3 w 4 for FEM and /3 « 8 for 
CGMD in rough agreement with the numerical solution; 
however, we should stress that this hypothesis has not 
been proved mathematically and the resonance structure 
of the scattering curve leads to large uncertainties in the 
/3 fit. If we fit to the tops of the peaks at long wave- 
lengths in scattering from an abruptly changing mesh, 
the uncertainty is reduced to /3 ~ 4 ± 0.2 and 8 ± 0.2 for 
FEM and CGMD, respectively. 

A series of bumps is visible in each of the curves in 
the transmissible region. Most of these bumps were not 



29 



FIG. 13: A plot ol the mesh used for the scattering calcula- 
tions. The atomic scale mesh, partially visible at the ends, 
extends infinitely to the right and left. The mesh is commen- 
surate with the underlying atomic lattice, and the largest cell 
size is Nmax = 20 atoms. 
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visible in plot of the reflection coefficients in Fig. 2 of 
Ref. with a linear scale. The log scale used in Fig. 
1121 brings out these features in regions with extremely 
low scattering. These bumps are scattering resonances, 
wavelengths at which the scattering cross section is in- 
creased because the frequency of the incoming wave is in 
resonance with an internal mode of the scattering region. 
Of course, they are much more peaked on a linear scale, 
where their width is an indication of the lifetime of the 
state. The curvature of the peaks in the log-linear plot 
is low, indicating short-lived resonances. The height of 
the peaks is an indication of the scattering strength. If a 
peak were high and narrow, it would indicate a strongly 
scattering localized mode, which would be pathological 
behavior in a concurrent multiscale simulation. In gen- 
eral, weak scattering with broad resonances (if any) is de- 
sirable. The wave reflections in CGMD a weaker and the 
resonances much less strong than in FEM, although the 
distributed mass FEM actually has a 10% higher thresh- 
old than CGMD because its frequencies are higher. 

It is also interesting to consider the transmission co- 
efficient, plotted in Fig. ^1 for CGMD and the two ver- 
sions of FEM. There is a simple relationship between the 
transmission and reflection coefficients, T = 1 — 7?., so in 
principle the calculation of one is equivalent to the other. 
However, because the log-linear plot brings out features 
near zero while suppressing features near unity, the two 
plots show different information. The transmission coef- 
ficient drops exponentially above the threshold, similar 
to quantum mechanical tunneling through a rectangular 
potential barrier or the transmission of evanescent waves 
in optics. As in those cases the transmission coefficient 
also decreases exponentially as the size of the scattering 
region is increased. One interesting feature of the trans- 
mission coefficient curves is the absence of resonances. 
The peaks are absent because the CG region lacks the 
degrees of freedom that would cause resonances at these 
high frequencies. 

We have calculated the scattering on many different 
coarse-grained regions. The general features of the reflec- 
tion coefficient curves remain as the mesh is varied, but 
the details change. One of the most pronounced changes 
happens if the mesh varies too abruptly. In this case, 
strong scattering resonances may occur near the thresh- 
old, even for CGMD, as shown in Fig. El Note the linear 
scale in this plot. The mesh used for the comparison be- 
tween FEM and CGMD, shown in Fig.^J was generated 
using a tanh function for the cell size to ensure smooth- 
ness. For comparison, we have plotted in Fig. El the 
reflection coefficient for a CG region of the same size but 



S 10 



10 



CGMD 
lump, mass FEM ■ 
20 dist. mass FEM ■ 







0.5 1 

Wave number k / k, 



1.5 







FIG. 14: A comparison of the transmission of elastic waves 
through a CG region in three cases: CGMD and two va- 
rieties of FEM. The dashed line marks the natural cutoff 
[ko = it/ {NjnaxO,)]- Notc that the bumps evident in the plot 
of the reflection coefficient are absent from the transmission 
coefficient. 
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FIG. 15: A comparison of the reflection of elastic waves from a 
CG region for a mesh with smoothly varying cell size and one 
with an abrupt change in cell size, both computed in CGMD. 
The dashed line marks the natural cutoff [fco — it/ (NmaxO.)]- 
Note that on the linear scale the resonances are not visible for 
the smooth mesh, but are quite pronounced for the abruptly 
changing mesh. 



consisting almost exclusively of cells of size N^ax = 20, 
and it is clear that the abrupt change in mesh size leads 
to much stronger resonances. The increased reflection of 
waves with k at resonance could lead to an unphysical 
size scale, and smooth meshes should be used to prevent 
this undesirable behavior. Apart from ensuring smooth- 
ness, we have not optimized the mesh, and it may be 
possible to reduce the scattering further still through a 
more optimized mesh. 

It should also be noted that the plots of wave reflec- 
tion on a log-linear scale provide a sensitive test of the 
numerical formulation. The CGMD data plotted in Figs. 
Elare not the same data plotted in Fig. 2 of Ref . ItR The 
data were noisy at the level of 10~^°, and we ultimately 
recalculated the stiffness matrix using Eq. (|38|l in order 
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to remove the noise. 

The reflection coefficient provides a measure to quan- 
tify the fD scattering properties of CGMD. What in- 
formation does it give about the reflection problem. The 
actual amount of reflection in any application is the prod- 
uct of the elastic wave power spectrum and the reflec- 
tion coefficient. Suppose the application of interest (e.g. 
crack propagation) can only tolerate scattering at 1% of 
the natural level. If 99% of the power is at wavelengths 
greater than Ac, then there will be an acceptable level of 
scattering if the mesh size is less than or equal to Ac. 

In applications like the crack propagation problem, it 
may be important to consider non-linear effects as well. 
In the anharmonic MD crystal, waves of sufficiently large 
amplitude will steepen into shock waves. The wave ve- 
locity increases with the pressure, so that a wavefront 
with a slow rise to a higher pressure will steepen into a 
step-like shock wave in which the abruptness of the rise 
is ultimately limited by dissipative processes at the front. 
As a result, compressive waves that are generated at the 
crack tip evolve into shock waves that have a strong im- 
pact on the crack if they return to it due to the boundary 
conditions. The reflection coefficient is a property of the 
system in the small amplitude, harmonic limit, and as 
such does not give any information about the behavior of 
shock waves. Shock waves, of course, have an abrupt rise 
and hence have power at short wavelengths localized at 
the wavefront. When a shock wave is incident on the CG 
mesh, the short wavelengths are reflected. Since the mesh 
spacing increases gradually, this reflection disperses the 
power at the front; i.e., the shortest waves are reflected 
flrst, then the next shortest and so on. The shock wave 
is dispersed and much of the power flows out to the CG 
mesh, so the reflected wave is a low amplitude wave that 
does not steepen into a shock wave. Thus, while some 
short wavelength components are reflected, they are not 
shock waves and do not appear to have an appreciable 
impact on the processes in the MD region. The majority 
of the power is carried out into the CG region, effectively 
delaying its return to the MD region by the transit time 
across the CG region. 

This dispersion and delay in wave reflection due to 
the CG region is the way CGMD and FEM/MD hy- 
brid motheods solve the reflection problem. Several other 
solutions have been proposed that make use of absorb- 
ing boundaries>22iSS2iS242£ Those techniques have much 
lower scattering of short wavelengths and hence a lower 
Kapitza resistance at zero temperature. They also in- 
volve considerable computer memory usage and consid- 
erable coding overhead. At this point it is clear that 
several approaches exist that solve the wave reflection 
problem in principle, but it is not yet clear which will of- 
fer the ease of use and scalability that will be demanded 
for widespread use in large-scale simulation. 



VIII. CONCLUSION 

Coarse-grained molecular dynamics provides a means 
of concurrently coupling MD with a coarsened descrip- 
tion of the mechanics similar to FEM. The practical for- 
mulation relies heavily on the properties of a crystal lat- 
tice, and it is therefore suited to solids. The formulation 
discussed here is based on a Hamiltonian, a conserved 
energy for the system, and is free from ghost forces. We 
have applied CGMD to three dimensional systems with 
interatomic potentials that are many-body in nature and 
extend well beyond nearest neighbors. In this article, we 
have elaborated on how CGMD is implemented to include 
anharmonic (non-linear) and flnitc temperature effects. 
It has been shown that these effects are described by an- 
alytic matrix formulas that may be precomputed prior to 
the simulation. The formalism is useful from the point 
of view of both providing the means to take the calcula- 
tions beyond the harmonic description and to be able to 
estimate the errors that are made in resorting to a har- 
monic model. In the process we have shown how CGMD 
is related to other concurrent multiscale methodologies 
such as the Quasicontinuum technique. 

We have also reported on several calculations done 
with CGMD in order to understand its properties. The 
elastic wave spectra for solid argon and tantalum have 
been computed in MD, CGMD and FEM for compari- 
son. In some cases it has been possible to derive largely 
analytic formulas for these spectra and to analyze their 
differences. In this investigation the MD spectrum has 
been taken as the ideal, and the CGMD and FEM models 
have each been formulated to give within its own frame- 
work the optimal agreement with MD. Several interesting 
results were found, including some reported briefly in pre- 
vious publications. First, both CGMD and FEM agree 
well with MD in the long wavelength limit, as expected. 
It was shown that CGMD provides a better description, 
in that the error is 0{k^), than FEM with errors of the 
order 0{k^). Second, throughout the Brillioun zone, the 
CGMD errors are smaller than those of FEM, especially 
at certain points on the zone boundary. Of course, the 
CGMD spectrum is exact in the limit of one atom per 
mesh cell. 

We have also used elastic wave spectra to examine 
the effect of the internal relaxation terms present in 
CGMD. This relaxation is the difference between forc- 
ing the atomic displacements to agree exactly with the 
interpolated displacement fleld and allowing fluctuations 
about the interpolated displacement fleld provided the 
interpolated fleld remains the best flt to the atomic dis- 
placements. We have shown how the relaxation effects 
can be eliminated to produce a rigid approximation to 
CGMD (R-CGMD), a formalism similar to the Quasi- 
continuum technique. We have compared the CGMD 
and R-CGMD elastic wave spectra and found that the 
two agree well in the long wavelength limit, but that 
CGMD provides a better description of elastic waves of 
short and moderate wavelengths. The internal relaxation 
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is able to give a better description of the energetics of 
waves changing rapidly on the mesh. 

Much remains to be done with CGMD, and we are ac- 
tively developing the model. Many questions of numer- 
ical efficiency still need to be addressed. A controlled 
means of bandwidth reduction for the CGMD stiffness 
matrix is needed. Also an efficient and consistent treat- 
ment of wave absorption is an open challenge. We have 
not discussed the implementation of CGMD on parallel 
platforms; the decomposition of the MD region into par- 
allel domains is straightforward, but the decomposition 
of the CG region is less obvious and is linked to the ques- 
tion of the stiffness bandwidth. Finally, the formulation 
of CGMD presented here provides the foundation for use 
in full-scale applications, a subject to which we plan to 
return soon. 
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APPENDIX A: ALTERNATE DERIVATION OF 
CG HAMILTONIAN 

In this Appendix we compute the CG Hamiltonian us- 
ing a more rigorous and straightforward, albeit laborious, 
approach. The choice to put this computation in an ap- 
pendix rather than the main text was made because it 
repeats a calculation done in the main text. It gives dif- 
ferent formulas for the mass and stiffness matrices, which 
are quite useful, and the positioning in an appendix is not 
intended to indicate that these derivations and formulas 
are less important than those in the main text. For this 
computation, as in Section llV Al we will concentrate on 
calculating the contribution of the potential to the par- 
tition function l|18(l : 

Zpot{uk;P) = JduAe-^P"-''^""- (Al) 

where once again we combine the atomic and spatial in- 
dices to form the SiVatom-dimensional configuration space 
for the atoms and the SA'nodo-dimensional space for the 
CG displacements. Notice that we have left A as a prod- 
uct of S functions, rather than using the Fourier repre- 



sentation as we did above: 

3W„odc 

A= n (A2) 

Suppose that the constraints were of a particularly simple 
form: 

3Ar„odc 

Asimpic^ -Q siu,-S,,u,) (A3) 

Then the evaluation of the integral would be easy. The 
first 3iVnodc atomic displacements would not be inte- 
grated, but rather just set to the corresponding Uj value. 
The remaining 3(iVatom — Nnodc) degrees of freedom 
would be integrated by completing the square and eval- 
uating the Gaussian integrals. The only problem is that 
the constraints are not of the simple form ljA3|l . 

We must introduce some mathematical formalism to 
transform the constraints into a form analogous to the 
simple kind ljA3|l . The approach we take in this deriva- 
tion is based on an explicit factorization of the 37Vatom- 
dimensional space into the direct product of the subspace 
spanned by the constraints and the orthogonal subspace. 
In order to accomplish this factorization we introduce 
two projection matrices in 3iVatom-dimensional space 

pir.'' ^ fjAfjxfkxy' fku (A4) 

= Nj^NjxNkxr' Nk, (A5) 
= N,J,, (A6) 

where repeated indices are summed, as usual, and the 
inverses are SiVnodo x SA^nodc matrix inverses. These are 
projection matrices in the sense that fju ~ fjn and 
P^fjv — 0. Since Nj^ is a linear combination of fku, 
it likewise holds that P^^^Nj^ = Nj^ and P^Nj^ = 0. 
The matrices are useful because P'^'^ projects onto the 
constrained subspace: 

pCG^^ = N,^f,,u, (A8) 
= N,^u, (A9) 

where we have used Eq. (jA6p . Thus, whenever P^^: g^^j-g 
on the configuration space its result depends only on the 
nodal displacements; it is completely independent of the 
unconstrained degrees of freedom. This is just what we 
need. 

Using the fact that the two projection matrices sum 
to the identity, S^i, = + P^^ , we can rewrite the 
potential part of the partition function 

Zpot = yrfuAe-^'5«-(p-+p,--)D,,(p-+p„-^)«„ 
= y"due-5'3("^-D^.''-+2".-D><^«M) X 

^g-^Pu,N,^Di_,^Nk„Uk (AlO) 
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where we have introduced the notation of the orthogonal 
and cross components of the dynamical matrix: 



— p^ n N- 



(All) 
(A12) 



Note that there is no D^^ defined or used in Eq. IjAlOp . 

We can now compute the integral (jA10|) with a slight 
trick. We can integrate over the constrained degrees of 
freedom 



practice, it is advantageous to choose e to be at the up- 
per end of the eigenspectrum of £>^^ so that D^^, is well 
conditioned. 

The same analysis can be carried out for the mass ma- 
trix. The result is 



M. 



jk 



(A20) 

— mNj^Ni^i^i for monatomic systems (A21) 
In this case, e is chosen to be the average mass. 



X e 2 



(A13) 



where the argument of the exponential is independent of 
the constrained subspace, so C is just a constant (inde- 
pendent of /3) and hence irrelevant. Now we can restore 
the constrained subspace in order to make the integral 
easier by inserting an integral that equals unity: 



\ -3JV„odo/2 

fl e-kP'-.^P^?^^ (AM) 

£[3 J 



where e is an arbitrary constant that we are free to de- 
termine below. This integral is unity because P^*^ is just 
the identity matrix on the constrained subspace. Insert- 
ing this expression into Eq. (|A13|I . we have 



-3Ar„odo/2 



(A15) 



where now we are left with an Gaussian integral over all 
space. Most importantly, the matrix in the Gaussian 



f) — ^ fP^'^ 



(A16) 



is non-singular provided e 7^ 0, as it must be since Zpot 
has been well defined at each step of the calculation. 

Now we complete the square to transform Eq. (jA15|) 
into a purely quadratic Gaussian integral using the shift 



(A17) 



we find 



Zpot — C 



^^-j3Ar„odc/2 
7 -\l/2 

fdetDj 



J) 



fA18) 



where 5N = A^atom — A^node- We find the important result 
that the CG stiffness matrix is given by 



(A19) 



and each matrix in this expression is well defined. In 
principle, this expression holds for any non-zero e; in 



APPENDIX B: EFFECTIVE POTENTIAL 

The CGMD formalism we have developed is an effec- 
tive theory in the sense that the short wavelength modes 
are integrated out to determine the effective interaction 
of the long wavelength modes. In field theory, an effec- 
tive potential is computed that is somewhat different in 
character. The typical approach would be to define the 
coarse-grained fields as an expectation value of the corre- 
sponding combination of microscopic degrees of freedom. 
Note that this differs from the approach we have taken in 
that we constrain the coarse-grained fields to equal that 
combination of microscopic degrees of freedom: they are 
identical and not an ensemble average. It is only the de- 
grees of freedom that are integrated out that are treated 
as an ensemble average. We have taken this approach be- 
cause at least in principle there can be bifurcation points 
in the trajectories of the coarse-grained degrees of free- 
dom that would be eliminated by defining them as en- 
semble averages. Nevertheless, the conventional effective 
potential approach has a certain theoretical elegance, and 
it could be useful in some contexts. We will therefore give 
a brief discussion of the CGMD effective potential and in- 
vestigate its usefulness. Hopefully, in the process we will 
eliminate any confusion that might arise between the two 
approaches. 

Again we start with the Helmholtz free energy in the 
presence of a source F(Jfe): 

Z(Jfe) = e-'^^'^-'"'^ = Jdudp exp (-/3i7^^^ - J,. • N^^u, 

(Bl) 

Just as in spin systems the magnetization is the deriva- 
tive of the Helmholtz free energy with respect to the ex- 
ternal field, the expectation value of the coarse-grained 
displacement is the derivative of this CGMD Helmholtz 
free energy with respect to the conjugate source: 

dF 

dJk 



= (Nkf^Uf,) = Ufc 



(B2) 



One could then go further and take the Legendre trans- 
form to find the Gibbs free energy G{uk) 



for which 



G = F + Uk-3k 



dG 
duk 



(B3) 



(B4) 
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and so is a minimum with respect to at equilibrium 
when Jfc — 0. 

The challenge is to derive an expression for G(ufc) given 
that it is expressed in terms of F{3k)- A calculation of 
Jfc(uj) is needed. Fortunately, the formalism is well de- 
veloped in field theory^^, and the result is that G{uk) 
is the generating function for one-point irreducible (IPI) 
Feynman diagrams. The IPI graphs are those that can- 



not be divided into two disconnected diagrams by break- 
ing a single internal line. This is a subset of the connected 
diagrams we have considered for CGMD (cf. Fig. [^J, so 
the effective potential approach does lead to a simplified 
formalism. We have not explored this model in depth. 
It might be interesting to do so, but we believe that the 
approach presented in the body of the text is more mean- 
ingful for coarse-grained solid mechanics problems. 
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